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FOREWORD 


This  report  is  one  of  two  computer  program  user's  manuals  prepared  by 
Aerotherm  Division  of  Acurex  Corporation  under  USAF  Contract  F04611-69-C-0081. 
Included  herein  is  Volume  I  of  the  manual  for  the  Aerotherm  Graphite  Surface 
Kinetics  (GASKET)  computer  code.  This  volume  describes  the  problems  solved  by 
the  code  and  presents  an  input  (card  format)  user's  guide  and  sample  problems. 
The  work  was  administered  under  the  direction  of  the  Air  Force  Rocket  Propul¬ 
sion  Laboratory  with  Mr.  Robert  J.  Schoner  as  the  Project  Officer. 

Mr.  John  W.  Schaefer  was  the  Program  Manager  and  Mr.  Mitchell  R.  Wool 
was  Program  Engineer.  The  GASKET  code  was  developed  by  Dr.  Kimble  J.  Clark. 

This  technical  report  has  been  reviewed  and  is  approved. 


A.  D.  Brown,  Jr.,  Lt.  Col.,  USAF 
Chief,  Technology  Division 
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ABSTRACT 


A  Fortran  IV  computer  code  is  described  which  computes  thermochemical 
ablation  rates  of  pyrolytic  graphite  surfaces,  as  a  function  of  surface  tem¬ 
perature  and  pressure,  assuming  all  heterogeneous  reactions  at  the  surface  are 
kinetically  controlled.  The  calculation  of  the  chemical  state  at  the  surface 
utilizes  a  film  coefficient  model  which  accounts  for  the  unequal  diffusion  of 
species.  All  homogeneous  reactions  between  gas-pha3e  species  adjacent  to  the 
surface  are  assumed  to  be  in  equilibrium.  The  specific  surface  reactions  which 
are  taken  to  be  kinetically  controlled  include  the  reactions  of  condensed-phase 
carbon  with  water  vapor,  carbon  dioxide,  and  hydrogen.  Several  options  are 
available  for  computation  of  the  boundary-layer  edge  state  required  in  the  dif¬ 
fusion  model,  including  isentropic  expansion  or  compression  and  normal  or  oblique 
shock  wave  calculation  procedures. 

This  computer  code  is  designated  as  the  Aerotherm  Graphite  Surface  Kinet¬ 
ics  (GASKET)  program  and  provides  surface  mass  balance  quantities  needed  for 
ablation  predictions  by  the  Aerotherm  Charring  Material  Ablation  (CMA)  code  or 
the  Aerotherm  Axisymmetric  Transient  Heating  and  Material  Ablation  (ASTHMA) 
code . 


Volume  I  of  this  report,  presented  herein,  contains  descriptions  of  the 
fundamental  physical  events  modeled,  the  mathematical  equations  solved,  the 
information  required  for  input,  and  the  results  output  by  the  computer  code. 

An  input  (card  format)  user's  guide  is  provided  along  with  sample  input  and 
output  listings  to  enable  an  unfamiliar  user  to  successfully  operate  the  code 
and  understand  the  results.  Volume  II  of  this  report  contains  supplemental 
information  on  the  specific  Fortran  IV  routines.  Included  are  program  listings, 
flow  charts,  and  definitions  of  Fortran  variables. 
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SECTION  1 
INTRODUCTION 

The  Aerotherm  Graphite  Surface  Kinetics  (GASKET)  computer  program 
is  described  in  this  user's  manual.  This  program  was  designed  specifically  for 
calculating  graphite  surface  thermochemical  ablation  rates  assuming  the  hetero¬ 
geneous  reactions  of  the  condensed-phase  surface  carbon  with  water  vapor,  carbon 
dioxide,  and  hydrogen  to  be  kinetically  controlled.  The  GASKET  code  was  developed 
from  the  Aerotherm  Equilibrium  Surface  Thermochemistry  (EST)  computer  program. 

As  described  in  Reference  1,  the  EST  program  is  extremely  versatile  (although 
it  excludes  kinetics)  and,  as  a  result,  the  preparation  of  input  data  decks 
is  rather  complex.  In  contrast,  the  GASKET  code  solves  a  much  more  restricted 
class  of  problems  and  accepts  a  much  more  simplified  form  of  input  data,  thereby 
reducing  considerably  the  amount  of  work  required  of  the  program  user. 

The  purpose  of  Volume  I  of  this  user's  manual  is  to  enable  an  un¬ 
familiar  user  to  effectively  utilize  the  Graphite  Surface  Kinetics  computer 
program.  To  this  end,  this  document  contains: 

•  a  brief  overview  of  the  purpose  and  capabilities  of  the 
GASKET  program  (This  Section) 

•  a  non-rigorous  discussion  of  the  theoretical  foundations 
of  the  calculation  methods  (Section  2) 

•  a  brief  description  of  the  numerical  solution  procedures 
and  their  FORTRAN  coding  (Section  3) 

•  a  specification  of  the  definitions  and  formats  of  the  input 
data  deck  (Section  4) 

•  several  examples  of  computer  program  input  decks  and  printed 
output  (Section  5) 

Volume  II  of  this  manual  contains  the  following  additional  program 
documentation : 

•  definitions  for  all  FORTRAN  variables  used 

•  listing  of  FORTRAN  IV  source  decks 

•  flow  charts  of  program  logic  for  each  FORTRAN  routine 
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Other  documents  of  potential  value  to  the  reader  requiring  a  more 
detailed  description  of  the  theoretical  fundamentals  of  the  GASKET  program  are 
References  2  and  3. 

:  l  i 

The  Graphite  Surface  Kinetics  program  computes  graphite  surface 
ablation  rates  under  the  assumption  that  the  condensed-phase  carbon  at  the 
surface  is  consumed  according  to  the  following  kinetically-controlled  reactions: 

C*  +  H20  ■*  CO  +  Hj  !  , 

C*  +  C02  J  2C0 
2C*  +  H2  ■*  C2H2 

where  *  denotes  the  condensed  phase  species. 

All  other  heterogeneous  reactions  between  the  surface  carbon  and  gas-phase  species 
in  the  boundary  layer  are  assumed  to  be  unimportant,  and  all  homogeneous  reac- 

I 

tions  involving  only  gas-phase  species  adjacent  to  the  surface  are  assumed  to 
be  in  equilibrium.  In  some  problems  the  elements  H  and/or  0  are  not  present  in 
the  edge  gas,  in  which  case  one  or  more  of  the  above  reactions  cannot  occur. 

GASKET  automatically  determines  which  of  the  three  reactions  are  possible  fn  a 
given  problem. 

The  GASKET  program  contains  a  number  of  options  which  can  be  used 
to  specify  the  boundary  layer  edge  state.  In  most  cases  the  program  of  interest 
is  the  chemical  response  of  graphite  nozzle  surfaces  or  graphite  ablative  models 

exposed  to  high-temperature  propellant  combustion  products.  The  sketch  below 

>  1  ! 

indicates  the  various  possible  edge  state  definitions: 

1.  From  a  specified  chamber  state,  a  one-dimensional,  real  gas, 
isentropic  expansion  through  the  nozzle  may  be  performed. 

The  edge  gas  state  at  amy  downstream  station,  eg.,  the  nozzle 
throat,  is  thus  determined l  . 

t 

2.  By  utilizing  the  solution  obtained  in  1.  and  specifying  a 
shock  angle  (0°  for  a  normal  shock,  as  is  the  case  here), 
the  static  conditions  behind  a  shock  may  be  calculated. 

3.  By  utilizing  the  solution  obtained  in  1.  and  specifying  the 
shock  angle,  the  isentropic  stagnation  conditions  behind  a 
shock  may  be  calculated. 

Inherent  in  the  edge  state  calculations  described  above  is  the  assumption  of 

,  i  i 

chemical  equilibrium.  Once  the  boundary  layer  edge  state  is  determined,  the 
GASKET  program  proceeds  to  compute  graphite  surface  ablatipn  rates  for  a  range 
of  surface  temperatures. 

i 


i 
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This  manual  attempts  to  acquaint  the  unfamiliar  user  with  the  GASKET 
program  by  providing  a  brief  summary  of  underlying  physical  principles,  a  speci¬ 
fication  of  the  input  mechanics,  and  some  sample  problems.  It  will  be  seen  that 
utilization  of  ,the  GASKET  code  is  quite  straightforward.  However,  instances 
will  occur  in  which  the  program  fails  to  solve  a  certain  problem  and  the  reason 
for  the  failure  will  not  be  apparent  from  the  manual.  These  cases  are  best  re¬ 
solved  by  direct  communication  with  the  program  or  manual  authors,  and  such  com¬ 
munication  is  encouraged. 
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SECTION  2 
THEORY 

This  section  briefly  summarizes  the  theory  upon  which  the  calcula¬ 
tion  methods  of  the  GASKET  program  are  based.  An  attempt  has  been  made  to 
present  these  fundamentals  in  a  fashion  most  appropriate  to  a  user  unfamiliar 
with  the  code.  For  this  reason,  and  in  the  interest  of  brevity,  explanations 
tend  to  be  nonvrigorous  and  perhaps  somewhat  heuristic.  No  attempt  is  made  in 
this  section  to  discuss  the  mathematical  methods  employed  to  solve  these  equa¬ 
tions.  The  solution  procedures  are  very  briefly  considered  in  Section  3  which 
is  concerned  primarily  with  the  program  coding. 

As  already  mentioned  in  Section  1,  the  GASKET  code  is  designed  speci¬ 
fically  for  calculating  thermochemical  ablation  rates  for  graphite  surfaces  in 
the  presence  of  certain  kinetically-controlled  chemical  reactions.  Such  calcula¬ 
tions  are  inherently  open  system  calculations,  since  the  relative  amounts  of 
chemical  elements  at  the  surface  depend  upon  various  mass  transfer  and  material 
degradation  rates.  In  general,  to  perform  these  ablation  calculations  the  thermo¬ 
chemical  state  of  the  gases  at  the  boundary  layer  edge  must  be  specified  so  that 
the  driving  potential  for  elemental  mass  diffusion  to  the  surface  can  be  estab¬ 
lished.  This  determination  of  the  edge  state  is  a  closed  system  calculation, 
since  the  relative  amounts  of  each  chemical  element  are  prespecif ied  (for  example, 
by  specifying  the  propellant  elemental  composition).  Furthermore,  in  performing 
such  closed  system  calculations  the  assumption  of  chemical  equilibrium  is  in¬ 
voked.  As  mentioned  in  Section  1,  the  GASKET  code  has  several  options  for  deter¬ 
mination  of  the  boundary  layer  edge  state,  including  real-gas  isentropic  expan¬ 
sion  or  compression  and  normal  or  oblique  shock  wave  calculations. 

It  is  apparent,  then,  that  the  GASKET  program  performs  calculations 
for  two  basic  types  of  systems:  1.)  equilibrium  closed  systems,  and  2.)  mixed 
equilibrium-nonequilibrium  open  systems.  The  former  are  discussed  in  Section  2.1 
and  latter  are  discussed  in  Section  2.2.  In  addition  to  calculating  the  chemical 
and  thermodynamic  states  of  the  two  basic  systems,  the  GASKET  program  also  cal¬ 
culates  and  prints  out  transport  properties  (e.g.,  viscosity,  thermal  conductivity, 
etc.)  appropriate  to  those  states.  The  methods  employed  by  the  program  to  cal¬ 
culate  these  transport  properties  are  briefly  discussed  in  Section  2.3. 
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2.1  CLOSED  SYSTEMS  IN  EQUILIBRIUM 

Closed  systems  have  been  defined  as  those  for  which  the  relative 
amounts  of  chemical  elements  are  prespecified.  In  the  following,  closed  systems 
are  (somewhat  artificially)  subgrouped  into  static  and  flowing  closed  systems. 

The  discussion  of  simple  static  closed  equilibrium  systems.  Section  2.1.1,  is 
particularly  significant  because  it  is  here  that  equilibrium  chemistry  relations, 
which  are  applicable  to  all  systems  subsequently  discussed,  are  first  introduced. 
Flowing  closed  systems  include  iser.tropic  expansions  or  compressions  (Section 
2.1.2)  and  flows  with  normal  or  oblique  shock  waves  (Section  2.1.3). 

2.1.1  Static  Closet.  Systems  in  Equilibrium 

Consider  K  chemical  elements,  N^»  introduced  into  a  previously  evac¬ 
uated  container.  In  general,  these  elements  v/ill  interact  to  form  a  number  of 
chemical  species*,  N.  (gas  phase)  and  (condensed  phases).  If  enough  time  has 
elapsed  so  that  thermodynamic  and  chemical  equilibrium  is  established,  the  thermo¬ 
dynamic  state  of  the  system,  including  the  relative  amounts  of  chemical  species 
present,  is  completely  determined  if  two  independent  thermodynamic  variables  are 
known**.  This  condition  may  be  stated  mathematically  by  examining  the  governing 
equations  for  such  a  system,  and  .  howing  that  the  number  of  independent  equa¬ 
tions  is  equal  to  the  number  of  unknown  quantities. 

Relations  expressing  the  formation  of  the  gaseous  chemical  species 
from  crhe  gaseous  chemical  elements  may  be  written  as  follows: 

K 

L  cki  \  in 

k=l 

Similarly,  formation  of  condensed  phase  species  from  the  gaseous  elements  is 
written: 


K 

°kJl 


k=l 


N,. 


N, 


(2) 


In  the  above,  represents  the  number  of  atoms  of  element  k  in  a  molecule  of 

species  i  (gas)  or  species  l  (condensed). 


"Chemical  species"  as  used  here  includes  molecular,  atomic,  ionic,  and  electron 
species . 

** 


Duhem's  Theorem. 
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If  the  gas  phase  species  are  assumed  to  individually  behave  as  ther¬ 
mally  perfect  gases,  then  the  equilibrium  relation  corresponding  to  reaction  (1) 
is 


k=l 


(T) 


or 

K 

ln  Pi  "  X  °ki  1X1  Pk  "  ln  KPi(T)  (3) 

k-1 

where  denotes  partial  pressures  and  (T)  is  the  equilibrium  constant  for 
the  formation  reaction  (1)  of  species  N^.  For  each  candidate  condensed  phase 
species 

K 

-  E  to  vt>  |4) 

k-l 


where 

=  indicates  the  existance  of  the  condensed  phase  species  in 
equili.'iri im  with  gas  phase  species,  and 

<  indicate  s  that  the  condensed  phase  species  will  not  be  present 
in  equilibrium. 


For  each  chemical  element  introduced  into  the  system,  the  conserva¬ 
tion  of  atoms  dictates  that  the  amount  of  any  element  k  ir.  the  gas  and  condensed 
phases  (regardless  of  molecular  configuration)  must  sum  to  the  total  amount  of 
element  k  in  the  system.  Mathematically,  this  may  be  written,  for  each  element 
k,  as 


Mass  fraction 
of  element  k 
input  to  the 
system 


IF  X  Vi 

i-1 


ft-1 


(5) 
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where  ^  is  a  composite  system  molecular  weight*  defined  by 

« - 1  t  % ♦ 1  •*  9^s°ofsrarm) 

i=l  £=1 

and  where  X^  is  a  mole  fraction  of  condensed  phase  species  i l  defined  as 

y  m  molecules  of  condensed  species  & 

*  *  total  gas  phase  molecules  i 

In  addition,  for  the  gas  phase  species,  there  exists  the  requirement  that  the 
partial  pressures  must  sum  to  the  total  system  pressure 

I 

^  Pi  =  P  (6) 

i=l 

Mixture  thermodynamic  properties,  such  as  specific  enthalpy,  are 
related  to  the  species  concentrations  by  equations  of  the  form 


h 


Pihi 


+ 


X  h 
l  l 


(7) 


Consider  now  the  numoer  of  independent  equations  for  the  system. 

The  number  of  gas  phase  equilibrium  relations  (3)  is  equal  to  the  number  of  gas 
phase  species  I  minus  the  number  of  elements  K  (because  equations  (3)  are  trival 
when  i=k) .  In  addition,  there  exists  a  relation  such  as  (4)  for  each  of  the  L 
candidate  condensed  phase  species  in  the  system.  Note  that  the  system  tempera¬ 
ture  is  contained  implicitly  in  equations  (3)  and  (4)  through  the  temperature 
dependence  of  the  equilibrium  constants.  There  are  K  conservation  of  elements 
equations  (5)  ,  one  for  each  atomic  element  introduced  into  the  system.  The 
requirement  that  the  partial  pressures  sum  to  the  system  pressure  (6)  contri¬ 
butes  one  additional  equation.  For  any  additional  thermodynamic  properties  of 
the  mixture  (enthalpy,  entropy,  etc),  there  exists  equations  such  as  (7). 

Consider  next  the  variables  appropriate  to  this  formulation  of  the 
problem.  The  relative  concentrations  of  the  I  species  in  the  gas  phase  are 


This  is  the  molecular  weight  appropriate  to  the  ideal  gas  equation  of  state  if 
condensed  phases  are  present. 
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given  by  the  P^'s  and  the  amounts  of  the  L  candidate  condensed  phase  species 
are  given  by  (most  or  all  of  which  may  be  zero) .  In  this  formulation,  the 
composite  system  molecular  weight,  ft  is  also  a  variable.  There  are  one  each 
of  the  mixture  thermodynamic  variables  T,  P,  h,  s,  etc.  The  number  of  variables 
and  available  independent  equations  may  be  summarized  as 


NO.  OF  SUCH 

EQUATION 

NO.  OF  SUCH 

VARIABLES 

VARIABLES 

NUMBER 

EQUATIONS 

Pi 

I 

(3) 

I  -  K 

L 

(4) 

L 

n 

1 

(5) 

K 

p 

1 

(6) 

1 

T 

1 

n 

of  the  type  (7) 

n 

total 

total 

variables 

I+L+n+3 

equations 

I+L+n+1 

Thus,  there  are  two  less  equations  than  there  are  variables,  and  so  if  two  in¬ 
dependent  variables  are  specified  (e.g.,  P  and  T)  in  addition  to  the  elemental 
composition,  then  closure  is  obtained  and  the  chemical  and  thermodynamic  state 
of  the  system  may,  in  principle,  be  determined. 


The  GASKET  program  performs  this  determination.  That  is,  to  deter¬ 
mine  the  equilibrium  thermodynamic  and  chemical  state  at  the  boundary  layer  edge, 
one  needs  only  to  furnish  the  GASKET  program  with  the  elemental  composition, 
the  candidate  gaseous  and  condensed  phase  species,*  and  two  thermodymanic  pro¬ 
perties.  One  of  these  properties  must  be  pressure,  and  the  otehr  may  be  either 
temperature,  enthalpy,  or  entropy.  Given  this  information,  the  GASKET  program 


In  order  to  perform  thermochemical  calculations,  certain  basic  thermodynamic 
data  (described  in  Section  4)  must  be  provided  for  each  candidate  species.  As 
one  of  the  options  of  the  GASKET  program,  the  user  must  decide  whether  or  not 
he  wishes  to  use  the  thermodynamic  data  built  into  the  program.  If  not,  he 
must  input  this  data  in  the  form  of  three-card  sets,  one  set  for  each  species. 
A  certain  amount  of  judgment  is  required  on  the  part  of  the  user  relative  to 
which  candidate  species  should  be  included  in  a  given  system.  Frequently  this 
judgment  is  avoided  by  simply  inputting  data  for  all  species  containing  com¬ 
binations  of  the  input  elements. 
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i  i 

,  .  i 

will  calculate  and  output  the  mold  fractions  of  each  candidate  species,  the 
temperature,  enthalpy,  entropy,  density,  effective  molecular  weight,  equili¬ 
brium  and  frozen  specific  heats,  the  isentropic  exponent,  and  a  few  other 
quantities  of  potential  interest.  Where  appropriate,  these  properties  are  out¬ 
put  for  the  gas  phase,  condensed  phase,  and  composite  system.  Also  output  are 
some  transport  properties  which  will  be  discussed  ih  Section  2.3.  The  exact 
convention  for  inputting  data  to  the  program  is  specified  in  Section  4,  and 

I 

some  example  proglems  of  this  type  are  presented  in  Section  5.  '  , 

2.1.2  Equilibrium  Isentropic  Flow 

The  relations  discussed  in  the  previous  section  for  static  equili¬ 
brium  may  be  applied  to  calculate  the  states  of  a  flowing  fluid  if  local  thermo¬ 
dynamic  and  chemical  equilibrium  is  assumed.  The  treatment  of  equilibrium  is¬ 
entropic  flows  is  discussed  in  this  section,  and  the  calculation  of  the  state 
downstream  of  a  shock  wave  is  discussed  in  Section  2.1.3. 

!  I 

For  adiabatic  flows,  the  conservation  of  energy  requires  that 

'  :  1 
h  +  I  u2  *  hQ  (8) 

i 

and  if  such  flows  are  in  equilibrium  (or  if  they  are  iully  frozen)  ,  they  are  , 
isentropic 

s  »  constant  (9) 

,  i 

I 

,  1  | 

The  application  of  these  relations  in  calculating  the  states  c}f  a  flowing  chemi¬ 
cally  reacting  g&s  may  be  demonstrated  by  considering  the  expansion  of  a  gas 

I 

from  a  plenum  chamber  through  a  nozzle.  In  the  plenum  chamber,  a  closed  system 
static  equilibrium  solution  (discussed  in  the  previous  section)  may  be  performed 
to  obtain  the  stagnation  conditions  (s,  hQ)  for  the  subsequent  isentropic  expan¬ 
sion.  Consider  next  the  state  at  some  point  in  the  nozzle.  Recall  that  in 
Section  2.1.1  it  was  shown  that  the  thermodynamic  and  chemical  state  of  a  closed 

I 

equilibrium  system  is  determined  by  two  independent  state  variables.  In  the 
present  case,  one  of  these  is  given  since  the  flow  is  isentropic  (9).  Thus, 
the  thermodynamic  and  chemical  state  of  the  fluid  may  be  completely  determined 
by  specifying  only  one  property,  conveniently  the  state  pressure,  P.  Once  the 
state  is  known,  the  local  velocity  may  be  calculated  from  (8). 


i 


I 


I 
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i 


i 


The  local  Mach  number  is  the  ratio  of  the  local  velocity  and  sound 
speed  where  the  local  sound  speed  is  given  by 


a  ■ 


j  [ 

For  simple  one-dimensional  internal  flow,  continuity  requires  that 


i 


(10) 


(11) 


In  (10):  and  (11),  o  is  the  total  system  density,  i.e.,  the  mass  of  gaseous  and 
condensed  phases  per  unit  volume.  When  performing  isentropic  flow  calculations, 

,  I 

the  GASKET  program  prints  out  (in  addition  to  the  local  thermodynamic  and  chemical 
state)  the  local  velocity  (8)  Mach  number  (10),  and  mass  flux  (11).  The  appropri¬ 
ate  input  for  isentropic  flow  calculations .  is  specified  in  Section  4,  and  some 
related t example  problems  are  presented  in  Section  5. 


2.1.3  Shock  Waves 

In  iadditio.1  to  equilibrium  isentropic  flows  (discussed  in  the  previ¬ 
ous  section) ,  the  GASKET  program  can  calculate  the  equilibrium  state  downstream 
of  a  normal  or  obiiqtie  shock  wave.  The  basis  for  these  calculations  is  discussed 
in  this  section.  ' 

i 

. 

To  illustrate,  consider  flow  through  an  oblique  shock  wave  with  up¬ 
stream  conditions  known  and  downstream  conditions  unknown. 


i 


i 


. 
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Equations  expressing  the  conservation  of  mass,  energy,  and  momentum  (indepen¬ 
dent  of  events  "in"  the  shock)  may  be  written  as  follows: 

Mass 

P2u2cos02  *  PjUj^osSj^  (12) 

Energy 

h2  +  7  u2  *  hl  +  7  U1  (13> 


Normal  Momentum 


P2  +  *>2u2cos 2 ®2  “  P1  +  P1^[co8iQ1 


y 


( 1 4  > 


Tangential  Momentum 


P2u|cos02sin02  ■  p^u^cosO^sinO^ 


(15) 


The  four  conservation  equations,  as  written  above,  involve  five  un¬ 
knowns:  02,  u2,  P2»  P2 »  h2.  However,  the  thermodynamic  variables  on  the  down¬ 
stream  side  of  the  shock  are  related  by  the  chemical  equilibrium  and  thermo¬ 
dynamic  state  relations  discussed  in  Section  2.1.1.  Recall  that  these  relations 
were  sufficient  to  completely  define  the  thermodynamic  and  chemical  state  of  a 
system  given  the  elemental  composition  and  any  two  independent  thermodynamic  state 
variables  (e.g.,  p2  as  a  function  of  P2  and  h2). 

Thus,  the  chemical  equilibrium  and  thermodynamic  state  relations, 
plus  the  four  conservation  equations  (12)  through  (15)  are  sufficient  to  determine 
the  conditions  downstream  of  a  shock  wave  if  all  upstream  conditions  are  known. 

The  GASKET  program  performs  this  determination  as  part  of  the  normal  or  oblique 
shock  option.  The  conditions  upstream  of  the  shock  may  be  specified  directly, 
or  they  may  be  carried  over  from  a  previous  isentropic  flow  solution  (Section 
2.1.2).  The  calculated  conditions  downstream  of  the  shock  include  all  thermo¬ 
dynamic  state  variables,  chemical  species  concentrations,  @2,  u,,  M2,  etc.  The 
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setup  o£  the  program  input  for  performing  a  shock  solution  is  discussed  in  Sec¬ 
tion  4  and  some  related  example  problems  are  presented  in  Section  5. 

2.2  MIXED  EQUILIBRIUM  -  NONEQUILIBRIUM  OPEN  SYSTEMS 

The  theory  upon  which  the  surface  ablation  calculations  are  based 
is  best  described  by  considering  a  sequence  of  physical  models  of  increasing  com¬ 
plexity.  Thus,  in  Section  2.2.1  both  simplified  and  more  refined  models  for  open 
systems  in  equilibrium  are  discussed.  Then,  in  Section  2.2.2,  further  refine¬ 
ments  to  the  open-system  model  to  account  for  chemical  kinetics  are  considered. 

2.2.1  Open  Systems  in  Equilibrium 

An  open  system  model  must  be  used  to  treat  thermochemical  ablation 
of  carbon  surfaces,  because  the  elemental  composition  in  a  control  volume  at¬ 
tached  to  the  surface  depends  upon  boundary  layer  transport  events  and  the  rate 
of  surface  material  degradation.  The  GASKET  program  uses  a  reasonably  sophis¬ 
ticated  approach  for  the  treatment  of  such  problems. 

As  will  be  discussed  in  more  detail  subsequently,  the  GASKET  pro¬ 
gram  relates  the  surface  state  (T  ,  h  ,  etc.)  to  the  ablation  rate.  The  surface 

w  w 

state  and  ablation  rate  are  also  related  by  surface  energy  balance  considerations 
which  depend  on  the  rate  of  heat  conduction  into  the  solid.  Thus,  the  particular 
ablation  rate,  surface  temperature,  etc.  is  calculated  by  combining  these  two 
relations:  surface  thermochemistry  +  in-depth  heat  conduction.  In  addition, 

the  boundary  layer  heat  and  mass  transfer  coefficients  must  be  specified.  Com¬ 
putationally,  the  surface  thermochemistry  relations  are  usually  generated  by  the 
GASKET  program  in  the  form  of  tables  (card  sets)  for  subsequent  input  to  tran¬ 
sient  heat  conduction  and  ablation  programs  such  as  the  CMA  (Charring  Material 
Ablation  -  Ref.  4)  program  or  the  ASTHMA  (Axisymmetric  Transient  Heating  and 
Material  Ablation  -  Ref.  5)  program.  The  information  flow  for  a  complete  tran¬ 
sient  ablation  prediction  may  be  illustrated  diagramatically  as: 


thermochemistry 
calculations 
(GASKET  program) 

C„  and  C„  -  B.L.  heat  and 
mass  transfer  coefficient 
distributions  &  histories 


in-depth  transient  heat  conduction  and 
ablation  calculations 
(CMA  or  ASTHMA  programs) 

surface  temperature  and  recession 
history,  in-depth  temperature  dis¬ 
tribution  history,  etc. 


surface  (and  perhaps  B.L. 
edge)  thermochemistry 
tables 


The  preparation  of  the  input  for  these  ablation  calculations  ia 
discussed  in  Section  4  for  the  GASKET  program  and  in  References  4  and  5  for  the 
CMA  and  ASTHMA  programs.  The  fundamentals  underlying  open  system  thermochemistry 
are  briefly  summarized  here.  Open  system  calculations,  subject  to  a  few  simplify 
ing  assumptions  to  more  clearly  illustrate  the  basic  theory,  are  discussed  first 
in  Section  2. 2. 1.1.  These  considerations  arc  extended  to  more  general  systems 
in  Section  2. 2. 1.2. 

2. 2. 1.1  Open  Systems  in  Equilibrium  -  Simplified  Case 

The  basic  theory  underlying  the  treatment  of  open  systems  may  best 
be  illustrated  by  examining  the  equations  expressing  the  conservation  of  chemical 
elements  and  energy  at  the  ablating  surface.  If  tho  boundary  layer  is  charac¬ 
terized  by  equal  diffusion  coefficients,  unity  Prandtl  number,  and  unity  Lewis 
number,  and  if  no  material  is  removed  from  the  surface  in  a  condensed  phase  (i.e. 
no  mechanical  erosion)  ,  then  these  equations  take  on  a  particularly  simple  form 
for  equilibrium  systems. 

Consider  first  the  fluxes  of  chemical  elements  (k)  entering  and 
leaving  a  control  surface  affixed  to  the  ablating  surface.  The  graphite  surface 
material  may  be  visualized  as  moving  into  this  surface  at  a  rate  s.  If  it  is 
assumed  that  no  material  is  being  remo”cd  in  a  condensed  phase,  then  the  surface 
and  the  fluxes  of  the  chemical  element  may  be  illustrated  as 


boundary 
layer  edge 


j 


w 


(cV>wKk 
.  w 


ablating 

surface 
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Terns  superscripted  by  a  tilde  (-)  represent  the  total  mass  fraction  or  flux  of 
element  k,  independent  of  molecular  configuration.  Thus 


I 


where  oertains  to  element  k,  i  pertains  to  species  i,  and  is  the  mass  frac¬ 
tion  of  element  k  in  species  i.  Fluxes  of  element  k  away  from  the  surface  con¬ 
sist  of  boundary  layer  diffusion  and  gross  motion  of  the  fluid  adjacent  to  the 
surface  due  to  the  injection  flux  ms.  Note  that  for  graphite  surfaces  the  quan¬ 
tity  K.  is  unity  for  k  «  C  and  zero  for  k  /  C  (where  C  represents  the  element 
*s 

carbon) . 

From  the  above  sketch,  it  is  seen  that  conservation  of  chemical 
elements  requires  that 


<»vluKk 


".'k 


s 


Summing  Equation  (16)  over  all  elements  k  yields  the  total  mass  continuity  equa¬ 
tion  (for  the  case  when  there  is  no  condensed  phase  material  removal) 


(PV»w  " 


m. 


(16) 


An  important  fundamental  of  the  present  mathematical  modeling  of  the  ablation 
process  is  the  expression  of  the  diffusive  heat  and  mass  fluxes  in  terms  of  a 
transfer  coefficient  formulation.  This  formulation  will  be  discussed  briefly  in 
the  following  paragraphs,  and  more  detailed  treatments  are  given  in  References 
6,  7,  and  8. 


2-12 


Heat  and  Man  Transfer  Coefficients  -  Simplified  Case 

For  low  speed  flow  of  an  incompressible,  constant  property,  non¬ 
reacting  fluid,  the  boundary  layer  energy  equation  may  be  written 
in  the  form 


u 


(17) 


and  solutions  to  this  equation  ))ave  historically  been  correlated  in 
terms  of  a  heat  transfer  coefficient  h  where 


l.(Te  -  T„> 


(18) 


For  high  speed  chemically  reacting  boundary  layers  (as  are  of  interest 
in  ablation  problems) ,  the  energy  equation  can  be  written  in  ' 
the  following  form  if  diffusion  coefficients  are  equal  and  if  t^he 
Prandtl  and  Lewis  numbers  are  unity  (e.g.,  Ref.  6) 


...  3H  3H  _  3_  /  3H\ 

pu  3j{  +  pv  3y  3y  3yj 


(19) 


where  H  is  the  total  (sensible  +  chemical  +  u*/2)  enthalpy.  By 
analogy  to  Equations  (17)  and  (18),  solutions  to  (19)  are  conveniently 
expressed  in  terms  of  a  dimensionless  heat  transfer  coefficient  CR 
where 

i 

«  “  peueCH(He  "  V  <20> 


Corresponding  formulations  for  non-unity  Prandtl  and  Lewis  number 
cases  will  be  considered  in  Section  2. 2. 1.2 

Consider  now  the  equation  for  mass  transfer  diffusion  in  the  boundary 
layer.  If  the  diffusion  coefficients  for  all  chemical  species  are 
equal  (Pick's  Law),  then  an  appropriate  summation  (Ref.  6)  of  the 
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boundary  layer  species  conservation  equation  yields  the  following 
equation  for  each  element  k 


(21) 


Again,  by  analogy,  solutions  to  this  equation  may  be  correlated  in 
terms  of  a  dimensionless  mass  transfer  coefficient  Cu  where 


'Ik  *  peueCM*\  '  Kk  1 
w  e  w 


(22) 


And  for  the  unity  Pradtl  and  Lewis  number  case  under  consideration, 
the  similarity  of  Equations  (19)  and  (21)  indicates  that  if  the 
corresponding  boundary  conditions  are  also  similar,  then  (?m  =  CH 

Utilizing  this  transfer  coefficient  formulation  (22)  for  the  diffusion 
flux  in  the  _lemental  balance  Equation  (16)  yields 


peueCM(Kk 


w 


\  1  *  lpv)wKk  -  "sV 
e  w  s 


(23) 


defining  a  dimensionless  ablation  rate  as 


m 


m 


(pv) 


B.  =  - 2L_ 

9  peueCM 


B '  s  - — r-  and 

c  peueCM 


B'  = 


w  _ 


m 


p  u  C..  p  u  C.. 

e  M  Fe  e  M 

(24) 


and  solving  (23)  for  the  total  mass  fraction  of  element  k  at  the  wall  yields  (for 
equal  diffusion  coefficients,  Pr  =  Le  =  1.0,  and  no  condensed  phase  material 
removal) 


w 


B '  K,.  +  K 


1  +  B' 


(25) 
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Given  the  relative  amounts  of  chemical  elements  specified  by  (25)  , 
the  chemical  and  thermodynamic  state  of  the  gases  adjacent  to  the  ablating  sur¬ 
face  may  be  calculated  from  equilibrium  relations  similar  to  those  discussed  in 
Section  2.1.1.  If  the  gases  are  in  equilibrium  with  the  ablating  surface 


in  P, 


k=l 


in  Kp£ (T) 


(26) 


if  i  represents  the  surface  species,  and 


-  Z  Sa  ln  pk  <  ln  vm 

k=l 


(27) 


for  all  other  candidate  condensed  phase  species.  For  a  graphite  surface,  however, 
it  is  assumed  that  _.ie  only  possible  surface  species  is  condensed-phase  carbon  and, 
further,  that  the  surface  is  kinetically  isolated  from  the  edge  gases.  Conse¬ 
quently,  GASKET  uses  rate  equations  in  place  of  equations  (26)  and  (27)  .  The 
equilibrium  relations  for  gas  phase  species  is  again  (see  Section  2.1.1  of  the 
form 


K 


in  P 


rls.i1"k  =  ,n  EPi(TI 


(28) 


k=l 


and  the  partial  pressures  must  obey  the  relation 


E'i- 

i=l 


(29) 


i 
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The  elemental  mass  fractions  adjacent  to  the  surface,  Kk  ,  of  (25)  are  related 
to  the  species  partial  pressures,  P^,  by  relations  such  as  (5)  (except  that  no 
condensed  phase  is  being  considered  here) 


c  P 
Sci  i 

i=l 


(30) 


Thus,  if  P  is  known  and  Tw  is  specified  (this  may  be  varied  parametrically  as 
will  be  discussed  subsequently),  and  if  B*  and  P^  are  unknowns,  then  the  number 
of  unknowns  and  equations  available  may  be  summarized  as 


UNKNOWNS 

NO.  OF  SUCH 
UNKNOWNS 

EQUATION 

NO. 

NO.  OF  SUCH 
EQUATIONS  AVAILABLE 

P. 

1 

I 

(28) 

I-K 

Kkw 

K 

(26) 

1 

V/( 

1 

(25) 

K 

B' 

1 

(29) 

1 

(30) 

K 

Total 

Unknowns 

I+K+2 

Total 

Equations 

I+K+2 

Thus,  closure  is  obtained  and,  in  principal,  the  carbon  ablation  rate 
(B')  and  molecular  composition  of  the  gases  adjacent  to  the  surface  may  be  det¬ 
ermined  if  P  and  Tw  are  specified.  From  the  pressure,  temperature,  and  chemical 
composition,  the  calculation  of  other  thermodynamic  properties  (enthalpy,  etc.) 
is  straightforward. 

The  GASKET  program  has  the  capability  to  determine  the  chemical  and 
thermodynamic  state  of  the  gases  adjacent  to  an  ablating  surface  in  a  fashion 
similar  to  that  discussed  here.  For  open  systems,  the  GASKET  output  consists 
of  card  sets  (in  addition  to  regular  printed  output)  referred  to  as  surface 
thermochemistry  tables.  Each  card  contains  information  relative  to  one  surface 
state  (i.e.,  one  P,  Tw,  B'...,  combination).  As  previously  discussed,  this 
surface  thermochemistry  data  alone  is  insufficient  to  solve  an  ablation  problem; 
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it  must  be  combined  with  in-depth  heat  conduction  and  surface  energy  balance  con¬ 
siderations.  In  order  to  more  clearly  demonstrate  this  interdependence,  the  re¬ 
lation  of  the  surface  thermochemistry  to  the  surface  energy  balance  is  discussed 
briefly  in  the  following  paragraphs. 

Relation  to  Surface  Energy  Balance  -  Consider  the  fluxes  of 
energy  entering  and  leaving  a  control  surface  affixed  to  the  ablat¬ 
ing  surface.  For  the  unity  Prandtl  number,  unity  Lewis  number,  no 
condensed  phase  removal  ablation  case  being  considered,  these  fluxes 
may  be  illustrated  as 


boundary 
layer  edge 


^conv  (pv)Hw 


ablating 

surface 


^cond 


**rad  in 


m_H 
s  s 


where : 

qconv  *8  the  ener*y  convected  to  the  surface  as  a  result  of  boundary  layer 

transport  events.  This  term  includes  effects  of  both  heat  conduc¬ 
tion  due  to  a  temperature  gradient  in  the  gas  adjacent  to  the  sur¬ 
face,  and  energy  transport  due  to  endothermic  or  exothermic  chemical 
reactions  at  the  surface.  For  the  equal  diffusion  coefficient,  Pr  *= 
Le  ■  1.0  case  under  consideration  here,  the  previously  discussed 
transfer  coefficient  form  of  this  term  is  (20) 

^conv  =  peueCH  (He  “  V  * 

(pv)wHw  is  the  energy  flux  leaving  the  surface  with  the  "blowing"  flux  (pv)w. 
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^rad  in 

ocV 

^cond 


is  the  incident  radiation  flux, 
is  the  surface  re-radiation  flux. 

is  the  heat  flux  conducted  into  the  solid  material,  g  .  ■  -  k 
„  _  _  k  3T\  COnd 

^cond  <Ty  I 

/w,  solid 

is  the  energy  flux  intering  the  surface  associated  with  the  surface 
ablation 


Conservation  of  energy  at  the  surface  requires  that 

'WVe  -  V  +  ”sHS  -  + 

ocT„ 


^rad  in 


^cond 


or,  utilizing  (24) 


“e  -  (X*B')HU  *  B'Hg  ♦ 


^cond 
peueCH  ‘w  "  peueCH 


oe  T  4 


(31) 


For  the  most  problems,  Hg  is  known,  PeueCH  is  independently  deter¬ 
mined  from  boundary  layer  calculations,  e  is  known,  qrad  is  known 
(frequently  zero) ,  and  Hg  is  a  known  function  of  T  .  Additionally, 

B'  and  Hw  are  related  to  Tw  through  the  surface  thermochemistry 
tables  generated  by  the  GASKET  program  as  previously  discussed. 

Based  on  cor siderations  thus  far,  all  quantities  in  (31)  except 
cjcond  are  determined  for  a  given  Tw«  The  in-depth  heat  con¬ 
duction  term,  qcomj»  depends  only  on  previous  heat  conduction  events 
and  the  current  Tw>  Thus,  it  is  seen  that  in  general  there  exists 
only  one  value  of  Tw  which  satisfies  (31)  .  This  Tw  is  usually 
determined  in  an  iterative  fashion  by  a  heat  conduction  and  abla¬ 
tion  program  (o.g..  Refs.  4  or  5)  which  utilizes  the  surface 
thermochemistry  tables  as  boundary  coi  tions.  Once  Tw  is  known 
for  a  given  time  in  a  transient  solution,  B',  surface  recession  rate, 
in-depth  temperature  distributions,  and  other  quantities  of  interest 
follow  directly. 
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The  above  has  been  a  brief  discussion  indicating  the  interdepen¬ 
dence  of  the  surface  thermochemistry  solution  and  the  in-depth  heat  conduction 
solution  for  a  somewhat  simplified  case.  It  should  be  pointed  out  that  simi¬ 
lar  considerations  apply  for  situations  when  unequal  species  diffusion  effects 
are  significant  and  when  Pr  f  Le  ft  1.0.  The  GASKET  program  is  able  to  treat 
all  of  these  effects  and  details  relative  to  the  treatment  of  some  of  these 
will  be  considered  in  Section  2.2.1. 2. 

2. 2. 1.2  Open  Systems  in  Equilibrium  -  Nonunity  Prandtl 

and  Lewis  Numbers,  Unequal  Diffusion  Coefficients 

The  discussion  of  the  previous  section  was  limited  to  open  systems 
with  unity  Prandtl  and  Lewis  numbers,  equal  species  diffusion  coefficients,  and 
no  removal  of  surface  material  in  the  condensed  phase.  These  simplifications 
were  made  in  order  to  render  the  basic  theory  easier  to  explain.  While  this 
simple  model  is  reasonably  accurate  for  many  ablation  situations,  these  assump¬ 
tions  are  inappropriate  for  others.  For  this  reason,  calculations  performed  by 
the  GASKET  program  (and  the  CMA  and  ASTHMA  programs)  are  not  restricted  to  any 
of  these  simplifications. 

The  equal  diffusion  coefficient,  Pr  =  Le  =  1.0  simplifications 
pertain  to  boundary  layer  mass,  energy,  and  momentum  transport  events,  and  the 
effects  of  the  relaxation  of  these  assumptions  on  the  problem  formulation  and 
solution  procedure  will  be  briefly  discussed  in  this  section. 

Pr  ?  Le  ft  1.0  -  For  nonunity  Prandtl  and  Lewis  numbers,  the  transfer 
coefficient  formulation  for  the  boundary  layer  energy  flux  is  not  as  straight¬ 
forward  as  that  discussed  in  the  previous  section.  This  is  because  the  boundary 
layer  energy  equation  is  no  longer  of  the  similar  form  and  thus  a  transfer  co¬ 
efficient  formulation  cannot  be  justified  purely  by  analogy.  A  detailed  dis¬ 
cussion  of  boundary  layer  transport  models  for  the  Pr  ft  Le  ft  1.0  case  is  beyond 
the  scope  of  this  manual  (see,  e.g.,  Refs.  7  and  8)  and  only  a  few  of  the  results 
as  they  relate  to  film  transfer  coefficient  formulations  will  be  discussed  here. 

When  the  Prandtl  number  is  not  unity,  the  viscous  dissipation  and 
heat  conduction  terms  in  the  boundary  layer  energy  equation  cannot  be  combined 
thus  rendering  the  equation  inhomogeneous.  Solutions  to  this  equation  indicate 
that  the  "driving  potential"  in  the  transfer  coefficient  expression  for  the 
surface  heat  flux  should  be  defined  in  terms  of  a  recovery  enthalpy  (e.g., 

Ref.  9)  in  place  of  the  actual  boundary  layer  edge  enthalpy.  Thus,  nonunity 
Prandtl  number  has  no  influence  on  GASKET  calculations  and  the  primary  practi- 
cal  implication  is  that  the  recovery  enthalpy,  rather  than  the  boundary  layer 
edge  enthalpy,  should  be  input  to  the  heat  conduction  and  ablation  solution. 
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When  the  Lewis  number  is  not  unity,  the  terms  in  the  boundary  layer 
energy  equation  representing  energy  transfer  by  heat  conduction  and  chemical 
species  diffusion  cannot  be  combined,  again  rendering  the  equation  inhomogeneous. 
The  energy  flux  to  the  surface  is  given  by 


4  «1  •  § 


h^°  pD 


(32) 


where  th<  first  term  characterizes  the  heat  conducted  to  the  surface  as  a  result 
of  the  temperature  gradient  in  the  gas  adjacent  to  the  surface,  and  the  second 
term  represents  the  effect  of  endothermic  and  exothermic  chemical  reactions  at 
the  surface.  The  appropriate  transfer  coefficient  form  of  (32)  is  not  firmly 
established  at  this  time.  However,  Reference  7  suggests  the  form 


P  uC 

e  e 


H(Hr  hwj frozen 
'  'edge  gas 


P  u  C 
e  e 


■  £  (\  -  \)‘> 


(33) 


and  this  is  the  form  employed  in  Aerotherm  ablation  programs.  In  (33),  the 

"driving  potential"  in  the  first  term  is  the  recovery  enthalpy  at  the  boundary 

layer  edge  minus  the  enthalpy  of  the  boundary  layer  edge  gases  frozen  at  the 

T 

edge  composition  and  at  the  surface  temperature,  and  h^w  represents  the  enthalpy 
of  chemical  species  i  with  respect  to  a  base  temperature  equal  to  the  surface 
temperature.  It  can  be  shown  (Ref.  7)  that  for  Le  =  1.0  and  CM  *  CH,  (33) 
collapses  to  (20)  as  expected.  However,  where  Le  /  1.0,  the  heat  and  mass  trans¬ 
fer  coefficients  are  generally  unequal  and  a  correlation  frequently  employed 
(Ral.  10)  is 


ffi 


Le 


2/j 


(34) 


It  is  apparent  from  (33)  that,  in  addition  to  surface  thermochemical  data,  bound¬ 
ary  layer  edge  thermochemical  data  are  also  required  when  Le  ?  1.0  in  order  to 
specify  the  quantity 
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“w  frozen 
edge  gas 


-I 

i-1 


Kiehi 


This  quantity  is  also  calculated  by  the  GASKET  program  for  subsequent  input  to 
heat  conduction  and  ablation  programs.  The  calculation  of  these  "frozen  edge 
tables"  is  usually  accomplished  by  performing  a  closed  system  equilibrium  cal¬ 
culation  at  the  boundary  layer  edge  temperature  and  elemental  composition  (as 
discussed  in  Section  2.1.1)  and  then  performing  closed  system  calculations 
frozen  at  this  composition  for  an  array  of  temperatures  spanning  the  expected 
wall  temperature  range.  The  program  input  for  frozen  composition  calculations 
is  discussed  in  Section  4.  Thus,  for  problems  with  CM  /  CH,  frozen  edge  tables, 
as  well  as  surface  thermochemistry  tables,  are  prepared  by  the  GASKET  program 
for  input  to  the  CMA  (Ref.  4),  ASTHMA  (Ref.  5),  or  related  heat  conduction  and 
ablation  programs. 

i 

Unequal  Diffusion  Coefficients  -  A  significant  simplification  of 
the  boundary  layer  energy  and  mass  diffusion  equations  results  if  all  binary 
diffusion  coefficients  for  a  given  species  i  are  assumed  to  be  equal  (e.g.,  Ref. 
9),  and  *1)  considerations  up  to  this  point  have  been  predicated  on  this  assump¬ 
tion.  However,  for  many  chemical  systems  of  interest  in  ablation  problems  (e.g., 
when  there  is  a  significant  difference  between  the  molecular  weights  of  the  majer 
species  present)  this  assumption  is  a  severe  compromise  with  reality.  Based 
upon  an  accurate  approximation  for  binary  diffusion  coefficients.  Reference  3 
presents  simplified  equations  tor  a  multicomponent  boundary  layer  with  unequal 
diffusion  coefficients  for  all  species.  The  application  of  similarity  arguments 
to  these  equations  suggests  the  following  form  for  the  transfer  coefficient 
formulation  for  the  diffusion  flux  of  element  k  at  the  surface 


(35) 


In  (35) ,  Z£  is,  in  effect,  a  weighted  average  of  the  mole  and  mass  fractions  of 
element  k.  The  Z£  are  defined  by 

*  i 


i 
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k  "  Z  V 

1-1 


n 


wYif  i“Y 
ZiKi 


1-1 


zV“Y 

ziKi 


(Y  =  2/3,  see 
Reference  2) 


VFi 

h  -  I - 

I]  vFi 

i-1 


where  the  F^  ere  diffusion  factors  defined  by  the  following  relation  for  the 
binary  diffusion  coefficients 


"Ij 


FiFi 


(36) 


where  D  is  a  constant  for  a  given  temperature  and  pressure  and  the  F^  depend 
weakly  on  temperature.  The  must  obey  (36)  in  order  for  the  boundary  layer 
species  diffusion  equations  to  reduce  to  a  form  from  which  (35)  can  be  inferred 
by  similarity  arguments.  Reference  3  demonstrates  that  the  binary  diffusion  co¬ 
efficients  for  a  varity  of  chemical  systems  are  accurately  correlated  by  (36). 
This  reference  also  shows  that  a  reasonably  good  correlation  equation  for  the 
F^  is 


where 


\ef 


23.4 


and 


e  «  0.431 


(37) 


I 


when  D  is  taken  as  the  self-diffusion  coefficient  of  Oj.  Additional  discussion 
relative  to  this  unequal  diffusion  coefficient  formulation  is  contained  in  Section 
2.3. 

•  | 

Consideration  of  unequal  diffusion  coefficients  affects  the  surface 
elemental  balance  relationships  which  are  needed  to  determine  the  equilibrium  state 


i 
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at  the  surface.  Substituting  (35)  into  (16)  yields 


peUeCM(Zku  “  Zk  1  +  (pv)wKk  "  msKk 
we  w  s 


and  the  "unknowns"  here  are  and  Z*  ,  each  of  which  may  be  expressed  in  terms 
of  the  species  partial  pressures 


E  Vi 


K-\  -T 


%  V 


E  Vi 


Lj 

9  i»l 


E  ‘Tel"'! 


ir*  1  =  1 
\  ¥ 


E  Sipi/FI  «Pi 

¥ - fV  E  c».ipi/FI 

E  VX 


where  is  the  mean  molecular  weight  of  the  gas  phase  and  F  is  a  mean  F?  defined 


E  Vi 

— 


Evi^l 

i-l 


Substituting  these  expressions  into  (38)  and  utilizing  (24)  yields  an  expression 
for  the  species  partial  pressures  at  the  surface  in  terms  of  quantities  at  the 
boundary  layer  edge  and  in  the  material 


*  Vn  P*  - 

B'  E  Vi*  * r  E  wl  ■  ¥  izk'  *  b’\. 

j 4  —  i  K  e 


2-23 


Note  that  (40)  reduces  to  (25)  when  the  diffusion  coefficients  are  equal. 

When  performing  unequal  diffusion  coefficient  open  system  calcula¬ 
tions,  the  GASKET  program  utilizes  (40)  rather  than  (25)  as  the  elemental  mass 
balance  equations.  Other  than  this,  the  solution  philosophy  is  essentially  as 
discussed  in  Section  2.2.1. 1.  The  diffusion  factors  utilized  in  the  GASKET 
program  may  be  calculated  in  three  ways,  at  the  user's  option 

a.  diffusion  factors  F^  may  be  input  individually  for  each 
species  i 

b.  diffusion  factors  may  be  calculated  according  to  (37)  with  the 
user  specifying  the  reference  molecular  weight, 7% re£ ,  and  the 
exponent  c 

c.  if  the  user  does  nothing  special,  the  program  will  calculate 
F^  according  to  (37),  with  2lref  K  23.4  and  e  ■  0.431. 

The  actual  program  input  for  these  alternatives  will  be  discussed  in  Section  4. 
It  should  also  be  pointed  out  that  the  diffusion  factors  have  an  effect  on  the 
other  transport  properties  calculated  and  printed  out  by  the  GASKET  program,  and 
these  will  be  briefly  discussed  in  Section  2.3. 

For  unequal  diffusion  coefficients,  the  transfer  coefficient  for¬ 
mulation  for  the  surface  energy  flux  has  the  form  (Ref.  7) 


WH«r 


* 

frozen  +  peueSi 
edge  gas  i*l 


(Zf 


w 


(41) 


Note  that  for  equal  diffusion  coefficients,  Z?  =  and  (41)  reduces  to  (33)  as 
expected.  Consistent  with  (41),  for  unequal  diffusion  coefficient  problems, 
the  surface  thermochemistry  tables  prepared  by  the  GASKET  program  contain  the 
quantity 


i-1 


w 


in  addition  to  the  quantities  previously  discussed.  Note  again,  that  for  equal 
diffusion  coefficients 
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V‘ 


(42) 


Similarly,  the  previously  diacussod  frozen  edge  table  must  contain  the  quantity 

Z*  h^w  «3) 

i-1  0 

in  addition  to 

hw  frozen 
edge  gao 

The  GASKET  program  input  arrangement  for  the  preparation  of  these  surface  and 
edge  thermochemistry  tables,  for  equal  or  unequal  diffusion  coefficients,  is 
specified  in  Section  4. 

2.2.2  Open  Systems  With  Kinotically  -  Controlled  Reactions 

To  calculate  the  equilibrium  state  of  a  chemical  system,  detailed 
atomic  data  for  each  possible  chemical  reaction  is  not  needed.  This  fact  per¬ 
mits  a  significant  simplification  in  the  problem  formulation,  and  the  theory 
briefly  outlined  in  Section  2.2.1  takes  advantage  of  these  simplifications. 

For  the  graphite  ablation  problem,  however,  it  has  been  shown  that  the  most  impor¬ 
tant  reactions  involving  the  condensed-phase  surface  carbon  and  gas-phase  species 
near  the  surface  do  not  satisfy  the  requirements  of  chemical  equilibrium  (see 
Ref.  11).  Thus,  an  accurate  treatment  of  the  thcrmochemical  ablation  of  graphite 
requires  an  open  system  model  which  accounts  for  both  equilibrium  and  nonequili¬ 
brium  chemical  reactions. 

A  general  solution  of  chemistry  problems  for  which  reaction  kinetics 
effects  are  important  is  potentially  difficult  for  at  least  two  reasons:  a)  there 
are  significant  computational  and  bookkeeping  problems  associated  with  the  analy¬ 
tical  treatment  of  mixed  equilibrium  and  nonequilibrium  systems,  and  b)  for 
many  systems  of  engineering  interest,  the  rate  controlling  reactions  are  not  well 
known  and/or  rate  constants  for  these  reactions  arc  unavailable.  The  GASKET 
program  effectively  surmounts  the  first  of  these  problems  as  discussed  in  Refer¬ 
ence  4  and  summarized  below.  With  regard  to  difficulties  falling  into  the  second 
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category,  experimental  and  analytical  atudiea  conducted  at  Aerotherm  have  resulted 
in  values  of  the  reaction  rate  constants  for  graphite  kinetics  which  are  the  most 
accurate  available  to  date  (References  11  and  12). 


In  the  remainder  of  this  section  the  GASKET  methodology  for  the 
treatment  of  chemical  kinetics  is  briefly  summarized.  In  Section  2. 2. 2.1  reac¬ 
tion  kinetics  in  general  are  discussed,  and  in  Section  2. 2. 2. 2  the  specific 
models  used  in  the  code  for  graphite  surface  kinetics  are  outlined. 


2. 2. 2.1  Reaction  Kinetics  in  General 

The  inclusion  of  kinetically  controlled  chemical  reactions  is  accomp¬ 
lished  by  removing  equilibrium  relations  such  as  equation  (3)  from  the  set  of 
equations  for  certain  specios  participating  in  reactions  that  are  to  be  kineti¬ 
cally  controlled.  These  equations  are  replaced  by  kinetic  rate  equations  for 
each  kinetically  controlled  reaction.  This  is  accomplished,  first,  by  identify¬ 
ing  the  primary  reactive  species  in  the  reactions  which  are  to  be  kinetically 
controlled  and,  second,  by  allowing  these  species  to  be  created  or  destroyed 
only  via  the  kinetic  rate  equations.  This  approach  is  performed  by  a  relabel¬ 
ing  of  species  to  be  considered  in  the  kinetically  controlled  reactions.  These 
species  are  called  pseudo-elements  since  they  behave  like  elements  except  that 
they  may  be  created  or  destroyed  at  rates  specified  by  the  reaction  rate  equa¬ 
tions. 

Computationally,  the  inclusion  of  kinetics  results  in  the  addition 
of  a  rate-of-creation  or  destruction  term  to  the  elemental  balance  equations  for 
these  pseudo-elements.  This  adds  additional  unknowns  to  the  system  equal  in 
number  to  the  number  of  species  whose  concentrations  are  kinetically  controlled, 
i.e.,  the  pseudo-elements.  The  relative  creation  and  destruction  rates  of  all 
pseudo-elements  in  a  given  reaction  are  related  by  stoichiometry,  however,  so  the 
number  of  additional  unknowns  remaining  is  equal  to  the  number  of  kinetically 
controlled  reactions.  The  reaction  rates,  from  which  the  pseudo-element  crea¬ 
tion  or  destruction  rates  derive,  are  given  by 


*Sn“  *Fi 


m 


TT/J"  -  k~  TT 


pm 


j 


pUjm 

Pj 


(44) 


for  each  kinetically  controlled  reaction  m  of  the  form 


1 


1 
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E^j  *  Eujr,Sj  (45’ 

j  j 

where  the  sums  and  products  are  over  the  species  and  pseudo-elements  N ^ ,  y  are 
the  stoichiometric  coefficients,  and  superscripts  R  and  P  denote  reactants  and 
products  respectively.  In  (44)  ,  Kpm  is  the  equilibrium  constant  for  reaction 
(45)  and  k^  is  the  forward  rate  constant.  In  the  present  formulation,  kFm 
is  represented  by  an  Arrhenius  type  function 

i 

$  I 

k  -  B  T  m  exp  (-E  /RT)  (46) 

Fm  m  ra 


where  Em  is  the  activation  energy,  is  the  temperature  exponent,  ahd  Bm  is 
factor  representing  a  multitude  of  phenomena.  These  quantities  are  generally 
based  on  experimental  data.  The  uncertainty  in,  or  unavailability  of,  these 
data  for  many  chemical  systems  of  interest  in  ablation  problems  frequently 
represent  a  significant  constraint  on  the  application  of  a  kinetics  model  to 
these  problems.  , 

When  certain  chemical  reactiorts  are  taken  to, be  kinetically  control¬ 
led,  the  arrangement  of  the  thermodynamic  data  for  all  candidate  species  must 
be  modified.  This  requirement  is  related  to  the  fact  that  all  species  involved 
in  kinetically-controlled  reactions  are  relabeled  and  designated  as  pseudo¬ 
elements.  Section  3.1.3  presents  a  more  complete  discussion  of  this  modeling 
subtlely  in  the  GASKET  program. 

2. 2. 2. 2  Kinetics  of  Hetergogeneous  Reactions  at  a  Graphite  Surface 

As  already  mentioned  in  Section  1,  in  the  graphite  ablation  problem 
the  condensed-phase  surface  carbon  is  kinetically  isolated  from  all  species  in 
the  boundary  layer  adjacent  to  the  surface.  That  is,  only  three1  species  are 
allowed  to  react  with  the  graphite  surface  material,  and  the  corresponding 
reactions  are  kinetically  controlled.  The  three  allowed  species  are  water  gas, 
carbon  dioxide,  and  hydrogen,  and  the  stoichiometric  equations  for  the  three 
reactions  are  as  follows: 

1)  C*  +  H20  ■*  CO  +  H2 

2)  C*  +  C02  J  2  CO 

3)  2C*  +  H2  l  C2H2  i 


l 
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The  detailed  development  of  the  graphite  kinetics  model  used  by  GASKET  can  be 
found  in  Reference  12.  Only  the  final  results  are  summarized  here. 

In  the  study  described  in  Reference  12,  it  was  found  that  the  effects 
of  reaction  inhibition  should  be  included  in  the  graphite  ablation  model.  That 
is,  experimental  evidence  indicates  that  the  kinetics  of  the  heterogeneous  sur¬ 
face  reactions  can  be  most  accurately  modeled  by  accounting  for  the  fact  that 
the  presence  of  certain  species  near  the  surface  tend  to  inhibit  the  forward 
tates  of  the  three  reactions.  The  method  of  Langmuir  (e.g..  Reference  13)  was 
used  in  the  development  of  this  heterogeneous  inhibition  kinetics  model.  The 
effects  of  reaction  inhibition  are  formally  treated  by  modifiying  the  pre¬ 
exponential  factor  in  Equation  (46): 

'  °'m  (i  *  t  ♦„  *,)(l  *  Z>J.  ' )  <4” 

j  j 

In  equation  (47)  B'm,  and  F^m  are  constants  for  the  Mth  reaction  and 
is  the  partial  pressure  of  the  j  species.  The  first  term  in  brackets  accounts 
for  the  inhibition  effects  due  to  any  species  which  is  also  involved  in  any  of 
the  three  reactions  above;  thus  the  summation  is  over  the  species  C* ,  I^O,  CO, 

H2*  Ct->2  an<*  C2H2*  The  second  term  in  brackets  accounts  for  the  inhibition 
effects  due  to  all  other  base  species  not  involved  in  the  three  kinetically 
controlled  reactions.  At  present,  the  quantity  Kjm  is  non-zero  only  for  the 
species  HC1  and  HF.  The  kinetics  model  comprised  of  Equations  (46)  and  (47)is 
referred  to  as  the  "Arrhenius/Langmiur "  model. 

The  inhibition  effects  are  assumed  to  be  the  same  for  layer- 
oriented,  edge-oriented,  and  bulk  pyrolytic  graphite.  The  matrix  of  values 
is  summarized  as  follows: 
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h2o 

CO 
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C2H2 

C* 

1 

1 

1 

1 

1 

0 

0 
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1 
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0 

3 

1 

1 

0 

1 

0 

0 
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The  Fjm  inhibition  effects  depend  upon  the  type  of  graphite  surface;  they  are 
the  same  for  bulk  and  edge-oriented  graphite  and  different  for  layer-oriented 
graphite.  The  matrix  of  K.  values  is  summarized  as  follows: 

1,2,  j  =  HCl ,  HF 

1,2,  3  =  all  species  other  than  HCl  and  HF 
3,  3  =  all  species 

The  other  constants  in  the  Arrhenius/Langmiur  model  which  are  built  into  GASKET 
are  summarized  as  follows: 


F .  =3  atm-1 ;  m  = 

]m 


m  = 


Surface 

lb  mole 

reaction/ft 1 sec 

cal/gm  mole 

unitless 

B'l 

B<2 

-  — , 

B*  3 

E1 

E2 

E3 

♦l 

*2 

*3 

Layer 

12.5 

12.5 

0.77 

46,000 

46,000 

55,500 

0 

0 

0 

Edge 

1.65x10- 

1.65x10- 

4 . 12x10  6 

65,500 

65,500 

129,500 

0 

0 

0 

Bulk 

1.98x10* 

1.98x10s  j 

1 . 94x10 7 

_ l 

65,500 

65,500 

L—  J 

129,500 

° 

0 

0 

These  latter  constants  were  also  determined  from  experimental  data,  as  described 
in  Reference  12. 

2.3  TRANSPORT  PROPERTIES 

In  addition  to  the  thermochemical  state  properties  discussed  in 
previous  sections,  the  GASKET  program  calculates  and  outputs  mixture  transport 
properties.  These  include  the  mixture  viscosity,  thermal  conductivity,  and 
species  diffusion  quantities.  The  Prandtl  and  Schmidt  numbers  are  also  output. 
These  transport  properties  are  calculated  from  expressions  which  derive  from 
simple  kinetic  theory  and  the  particular  multicomponent  diffusion  representation 
previously  discussed  in  Section  2. 2. 1.2.  The  development  of  these  expressions 
is  discussed  in  detail  in  Reference  3.  A  brief  summary  of  this  development,  and 
the  resulting  expressions,  are  presented  in  this  section. 

Diffusion  Coefficients  -  In  Section  2. 2. 1.2  a  bifurcation  approxi¬ 
mation  for  binary  diffusion  coefficients  was  mentioned  which  characterized  multi- 
component  diffusion  phenomena  with  reasonable  accuracy  without  unduly  complica¬ 
ting  the  system  of  equations  to  be  solved.  This  simplification  is  achieved 
through  a  correlation  for  binary  diffusion  coefficients  of  the  form 
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FiFj 


(48) 


where  5  is  a  reference  diffusion  coefficient  and  the  are  diffusion  factors. 
These  quantities  are  discussed  in  detail  in  Reference  3.  The  incorporation  of 
(48)  in  the  Stefan-Maxwell  relation  for  mass  diffusion  fluxes  indicates  that 
the  diffusion  flux  of  species  i  may  be  written  in  terms  of  only  properties  of 
species  i  and  global  system  properties.  Subject  to  a  few  simplifying  assump¬ 
tions  (Ref.  3),  this  expression  for  may  be  written 

pPp2  3zi  (49) 

3i  “  ft p  3y 


where 


*i  **  ixi/Pi>‘2 


(50) 


“1  -  Z“jF) 

j 

"2  * 

3 


(51) 


(52) 


The  accuracy  of  this  formulation  is  examined  in  Reference  3  for  a  variety  of 
chemical  systems.  It  is  shewn  that  the  2*^  calculated  by  (36)  represent  a  very 
substantial  improvement  over  equal  diffusion  coefficients  when  compared  to  exact 
values  calculated  directly  from  kinetic  theory.  The  calculation  of  the  mixture 
viscosity  and  thermal  conductivity  is  based  on  the  diffusion  factors  given  by 
(48),  and  these  will  be  discussed  in  the  following  paragraphs. 

Mixture  Viscosity  -  The  expression  employed  by  the  GASKET  program 
to  calculate  the  mixture  viscosity  derives  from  rigorous  first  order  kinetic 
theory  (Ref.  14),  subject  to  a  few  simplifying  assumptions,  as  discussed  in  Refer¬ 
ence  3. 
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V1*  385  P*L  - 

1  3=1  13 


,•  ^  a  a 

3=1 

j^i 


(53) 


where  ^  is  the  viscosity  of  the  pure  species  i.  The  ^  may  be  expressed  in 

terms  of  the  self  diffusion  coefficients  £  .  . 

n 


pi 


6 pi*ii 


(54) 


where  A*^  is  a  ratio  of  collision  integrals  based  on  a  Lennard-Jones  intermole- 
cular  potential.  Substituting  (48)  and  (54)  into  (53)  results  in  the  following 
expression  for  the  viscosity  of  the  multicomponent  mixture. 


mix 


fiD 
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1  i-1 


x^i 


1.385  + 


XiFi  /6Aii 

“l-5- 


-  1.385 


i 


(55) 


and  this  is  the  expression  utilized  to  calculate  the  mixture  viscosity  output 
by  the  GASKET  program. 


Mixture  Thermal  Conductivity  -  The  thermal  conductivity  in  a  polya¬ 
tomic  gas  mixture  may  be  represented  by  (Ref.  14) 

kmix  ~  kmono-mix  +  kint 

where  ^mono_mix  Is  the  thermal  conductivity  in  a  mixture  computed  neglecting  all 
internal  degrees  of  freedom  and  kinfc  is  the  contribution  to  the  thermal  con¬ 
ductivity  of  the  mixture  due  to  the  internal  degrees  of  freedom  of  the  molecules. 
A  simplified  expression  for  the  mono-mixture  thermal  conductivity  can  be  derived 
in  a  manner  similar  to  the  procedure  previously  discussed  for  the  mixture  vis¬ 
cosity.  This  simplified  expression  is  (from  Ref.  3) 
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where  mono  is  the  thermal  conductivity  of  the  pure  species  i  neglecting  all 
internal  degrees  of  freedom  of  the  molecule.  The  mono  may  be  expressed  in 
terms  of  the  y,  as  per 
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k 


i 


mono 


15  JL 
T  57^ 


(58) 


The  contribution  to  the  thermal  conductivity  from  the  internal  degrees  of  free¬ 
dom  may  be  expressed  as  (from  Ref.  3) 


^int 


E 


(59) 


By  combining  (48)  with  (56)  through  (59),  the  mixture  thermal  conductivity  may 
be  written  as 


where  y^  and  y2  are  given  by  (51)  and  (52)  respectively,  and 
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Thus,  (60)  is  the  expression  utilized  to  calculate  the  mixture  thermal  conducti¬ 
vity  output  by  the  GASKET  program. 


Also  calculated  and  output  by  the  GASKET  program  are  the  Prandtl 
and  Schmidt  numbers  which  are  defined  here  as 


Pr 
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Jc  p-frozen 


(63) 
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The  transport  properties  calculated  by  the  GASKET  program  are  all 
based  on  the  bifurcation  approximation  for  the  binary  diffusion  coefficients 
expressed  in  (48), .  This  is  so  even  for  closed  system  calculations  (in  which 
case  diffusion  phenomena  need  not  be  considered  to  calculate  the  chemical  and 
thermodynamic  state  of  the  system)  and  for  open  system  calculations  for  which 
equal  diffusion  coefficients  are  assumed  (Section  2. 2. 1.1).  From  the  equations 
presented,  it  may  be  observed  that  the  properties  calculated  are  highly  dependent 
on  the  diffusion  factors,  F^.  Three  alternate  methods  for  prescribing  the  F^ 
were  discussed  Section  2. 2. 1.2.  The  use  of  the  diffusion  factor  correlation  (48) 
with  resident  values  °^^cef  an< *  e  (which  were  derived  primarily  from  considera¬ 
tion  of  species  diffusion  coefficients)  should  result  in  reasonably  accurate 
values  of  other  transport  properties.  Alternately,  the  correlation  (48)  may  be 
used  with  values  of  and  e  derived  by  correlating  available  transport  prop¬ 

erty  data  for  the  particular  system  of  interest.  If  transport  properties  of 
maximum  accuracy  are  to  be  calculated,  then  the  diffusion  factors  should  be  in¬ 
put  individually  for  each  species.  These  data  may  be  obtained  from  tabulations 
such  as  Reference  14. 
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SECTION  3 
CODING 

This  section  presents  some  details  relative  to  the  FORTRAN  coding 
of  the  GASKET  program.  A  brief  description  of  the  numerical  solution  procedure, 
which  includes  a  discussion  of  the  "base  species"  concept,  is  presented  first 
in  Section  3.1.  Each  subroutine  making  up  the  GASKET  program  is  identified  and 
briefly  discussed  in  Section  3.2,  and  definitions  of  several  output  quantities 
are  provided  in  Section  3.3.  Some  miscellaneous  program  details  are  documented 
in  Section  3.4.  These  include  storage  requirements,  tape  requirements,  operator 
controls,  dumps  and  diagnostics,  and  run  time. 

3.1  COMPUTATIONAL  PROCEDURE 

The  computational  procedures  employed  in  the  GASKET  program  to  solve 
the  equations  set  forth  in  Section  2  are  briefly  discussed  here.  Considerably 
greater  detail  is  presented  in  Reference  2  and  in  particular  Table  I  of  that 
reference. 


3.1.1  Basic  Solution  Technique 

The  basic  solution  technique  may  be  illustrated  by  considering,  for 
example,  an  open  system  with  unequal  diffusion  coefficients,  no  condensed  phase 
material  removal,  and  no  kinetics  (i.e.,  as  discussed  in  Section  2. 2. 1.2).  For 
this  system,  the  basic  equations  defining  the  problem  are  the  elemental  conserva¬ 
tion  equations  (40) ,  the  total  pressure  equation  (29) ,  the  reaction  equilibrium 
equatiors  (28),  and  one  heterogeneous  vapor  pressure  relation  (26).  The  table 
of  Section  2. 2. 1.1  shows  that  there  are  many  knowns  as  unknowns  in  these  equa¬ 
tions  so  closure  is  obtained. 


Summarizing  these  equations, 
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I  ^  P19L 
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(67) 


-  E  V"  h  -  tn  VV  -  ®. 
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(68) 


The  number  of  unknowns  could  immediately  be  reduced  by  I-K  through 
the  direct  substitution  of  (67),  as  solved  for  P^  into  the  other  relations.  If 
proves  advantageous,  however,  to  continue  to  treat  the  full  set  of  equations, 
and  to  subsequently  utilize  this  substitution  during  the  iterative  convergence 
procedure.  The  solution  of  these  simultaneous  nonlinear  algebraic  equations 
is  based  on  Newton-Raphson  iteration.  Since  this  procedure  is  accelerated  by 
casting  the  equations  into  a  more  linear  form  (via  transformations,  substitu¬ 
tions,  etc.)  it  is  well  to  examine  the  set  of  equations  above.  With  the 
boundary  layer  edge,  char  and  pyrolysis  gas  composition  given  as  well  as  the 
B' ,  (65)  and  (66)  are  linear  relations  between  the  P|  and  171  providing  that 
F  is  reasonably  constant.  In  contrast,  (67)  and  (68)  are  linear  relations  bet¬ 


ween  the  In 
linear  in  1/T. 


tn  P.  and  Jin  K  . 

K  pi 


the  latter  variable  being  approximately 


The  GASKET  program  takes  advantage  of  this  situation  by  treating 
those  species  which  are  significant  in  the  mass  and  pressure  balances  (65) 
and  (66)  in  terms  of  P^  and  the  lpss  significant  species  in  terms  of  their 
*n  Pr  , 

The  Newton-Raphson  procedure,  as  applied  by  the  GASKET  program, 
can  be  summarized  by  considering  a  set  of  equations  of  the  general  form 


f j  (x1#x21 


• • i , . . . )  —  0 


i 


At  any  point  in  the  solution  procedure  there  exists  a  set  of  estimates,  xj, 
for  all  the  variables  which  will  in  general  not  satisfy  all  of  the  relations 
and  will  lead  to  a  non-zero  value  of  the  fj,  namely,  e  j .  The  Newton-Raphson 
method  proceeds  to  "drive”  these  errors  toward  zero  by  evaluating  the  change  in 
each  unknown  variable,  Ax^,  which  would  reduce  all  the  errors  to  zero  if  the 
functions,  f y  were  linear.  The  linear  approximation  is  based  on  the  current 
values  of  the  unknown  variables  and  the  corresponding  array  of  values  of  the 
partial  derivatives  3fj/3x^.  Thus 


(69) 


which  is  locally  correct  and  is  integrated  to 


(Af  j)  * 


ip  the  linear  approximation.  The  solution  of  (69)  is 


i 


dx. 


(70) 


where  the  array  of  partial  derivatives  appearing  in  (70)  is  simply  the  matrix 
inverse  of  the  array  in  (69).  >  In  the  GASKET  program  the  formulation  of  the 
partial  derivatives  uses  the  variables,  in  P^,*  in  T  and  lnU\,  and  (70)  yields, 
for  example,  , 


d(in  P^) 


?  in 


(71) 


which  if  taken  as  linear  all  the  way  to  solution  yields 
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since  the  desired  change  in  the  functions  is  simply  the  negative  of  the  error. 
An  equally  exact  relation  obtained  from  (71)  is 
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E3  in 
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df . 


which  if  taken  as  linear  all  the  way  to  solution  yields 


Pt*  ~  Pi  “  AP 
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(73) 


The  GASKET  program  uses  (73)  for  all  species  significant  in  mass  balances  and 
(72)  for  the  others. 

3.1.2  Restriction  on  Corrections 

* 

The  set  of  correction  (Ax^)  can  be  thought  of  as  a  vector  in  the 
space  of  the  independent  variables  which  is  added  to  the  current  vector  approxi¬ 
mation  xj  to  yield  a  new  estimate  x£*.  Experience  has  shown  that  it  is  fre¬ 
quently  unwise  to  proceed  along  this  correction  vector  the  full  amount  indicated 
by  (72)  or  (73).  Rather,  it  is  better  to  proceed  a  limited  way,  although  pre¬ 
serving  the  same  direction.  At  other  times,  it  is  expedient  to  depart  from 
this  vector,  and  seek  another  based  on  freezing  the  value  of  some  variable  and 
eliminating  a  corresponding  equation. 


The  choice  of  In  P.  permits  a  matrix  reduction  by  the  use  of  the  simple  alge¬ 
braic  substitution,  previously  mentioned  after  (68) ,  prior  to  matrix  inversion. 
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The  scaling  of  the  correction  vector  is  such  as  to  limit  changes 
in  the  partial  pressures  of  major  species  to  increases  of  one  order  of  magnitude 
and  decreases  of  three  orders  of  magnitude,  and  changes  of  temperature  to  approxi¬ 
mately  20  percent. 

Molecular  weight,  temperature  and  condensed  phase  concentration 
corrections  are  frozen  and  a  new  correction  vector  generated  if  the  initial  set 
of  corrections  indicate  excessive  temperature  or  molecular  weight  excursions, a 
contradictory  temperature  change,  or  negative  corrections  on  newly  introduced 
condensed  species. 

The  formulation  of  these  and  other  scaling  and  freezing  criteria 
is  an  essential  feature  of  the  GASKET  program.  Because  of  these  features,  con¬ 
vergence  is  virtually  assured  for  well  formulated,  physically  unique  problems. 

3.1.3  Base  Species 

The  discussion  of  Section  2  described  the  equilibrium  ruju  ;un 
equations  as  equations  giving  the  formation  of  a  species  from  atomic  elements. 

Thus  the  reactants  are  elements  and  the  products  are  usually  molecules.  This 
scheme  has  the  advantage  of  formal  simplicity,  since  the  stoichiometric  coef¬ 
ficients  needed  in  the  equlibrium  equations  are  given  directly  by  the  product 
species  chemical  formula.  This  scheme  can  have  computational  disadvantages, 
however,  since  the  atomic  elements  are  frequently  not  present  to  any  great 
extent  in  the  equilibrium  system.  This  in  itself  results  in  no  disadvantage. 

If,  however,  a  molecule  (e.g.,  CO)  dominates  more  than  one  mass  balance  (e.g., 

C  and  0)  loss  of  significant  figures  can  slow  or  defeat  convergence. 

It  is  more  desirable  to  write  the  equilibrium  reactions  (i)  and 
(2)  as  well  as  the  mass  balances  in  terms  of  reactant  species  which  are  actually 
present  in  appreciable  amounts.  These  species  are  termed  "base  species"  (from 
Ref.  15)  since  all  other  species  are  taken  to  be  formed  from  them. 

The  GASKET  program  selects  the  base  species  from  the  candidate 
species  thermochemical  data.  The  program  automatically  selects  as  base  species 
the  first  set  of  species  satisfying  the  requirement  that  (1)  all  other  species 
may  be  formed  from  this  base  species  set  and  (2)  that  no  balanced  reaction 
can  be  written  involving  only  base  species.  One  base  species  may  be  considered 
to  represent  each  element.  Thus,  the  base  species  are  established  by  the  order 
of  the  candidate  species  thermochemical  data.  The  calculation  of  the  stoichio¬ 
metric  coefficients  and  equilibrium  constants  appropriate  to  any  set  of  base 
species  is  handled  automatically  by  the  program. 
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As  mentioned  in  Section  2. 2. 2.1,  when  chemical  kinetics  are  con¬ 
sidered  the  arrangement  of  tho  candidate  species  thermochemical  data  must  be 
modified,  because  all  species  involved  in  kinctically  controlled  reactions  are 
identified  as  pseudo-elements.  The  treatment  of  chemical  kinotics  in  CASKET 
requires  tha.  each  of  the  species  involved  in  kinotically  controlled  reactions 
also  be  identified  as  a  base  specie.  Therefore,  tho  thormochemical  data  must 
be  arranged  so  that  the  species  C",  ll20,  CO,  «2,  COj  ,  and  Cjllj  appear  at  tho 
top  of  the  list  and,  thus,  can  be  selected  as  base  species  when  required.  The 
thervochemical  data  built  into  GASKET  is  arranged  in  this  required  ordor. 

Another  subtle  aspect  of  the  base  specie  selection  procedure  occurs 
when  chemical  kinetics  is  considered.  Tor  certain  edgo-gas  compositions  the 
number  of  elements  involved  is  fewer  in  number  than  tr.e  number  of  species  pre¬ 
sent  in  the  three  kinetically-controlled  reactions.  When  this  occurs,  additional 
"elements"  must  be  created  in  order  that  the  number  of  corresponding  base  species 
is  large  enough  to  include  all  species  in  the  reactions.  An  example  will  serve 
to  illustrate.  Suppose  that  the  edge  gas  contains  II,  C,  and  O.  If  only  these 
elements  are  considered,  then  only  three  base  species  arc  possible  and,  there¬ 
fore,  only  three  of  the  species  C*  ,  11^0,  CO,  ,  CO^ ,  and  C2H2  appearing  in  the 
kinetically-controlled  reactions  can  possibly  be  base  species.  This  situation 
is  remedied  by  taking  two  steps.  First,  CASKET  defines  two  new  "elements  C*  and 
C02 ,  identified  by  fictitious  atomic  numbers,  whose  corresponding  base  species 
are  C*  and  COj,  respectively.  Those  new  "elements"  arc  added  to  the  sot 
comprised  of  H.  C,  and  0.  Then,  GASKET  reinterprets  the  kinctically  controlled 
reaction  3)  as 

2C*  ♦  H2  -  C2  h2  -  2C0  ♦  3H2  -  2H20 

That  is,  the  products  of  a  kinctically  controlled  reaction  can  be  replaced  by 
their  equilibrium  base  spccicB  in  stoichiometrically  correct  proportions.  Thus, 
with  the  "elemental"  set  comprised  of  H,  C,  0,  C02  ,  and  C*,  five  base  species 
are  possible.  And,  with  the  thenuochcmical  data  in  the  correct  order,  the  fiv^ 
base  upee'es  selected  will  be  H2,  CO,  HjO,  CO  ^ ,  and  C* . 

Finally,  one  rubtletv  remains  to  be  considered.  In  the  previous 
paragraph  it  was  r.  jntioned  that  the  "element"  C*  was  defined  and  assigned  a 
fictitious  atomic  number.  It  follows  that  C*  can  communicate  with  the  remainder 
of  the  chemical  system  only  through  the  three  kinetically-controlled  reactions. 
With  this  in  mind,  consider  the  second  fictitious  "element"  COj.  In  reality, 
it  is  possible  that  C02  can  be  in  equilibrium  with  many  other  gas  phase  species 
while  simultaneously  reacting  with  the  graphite  surface.  The  only  way  this  can 
be  accomplished  is  to  introduce  explicitly  a  kinetically-controlled  reaction 
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involving  CC^  and  other  gas  phase  species  and  then  allow  the  rate  of  this  reac¬ 
tion  to  approach  infinity,  i.o.,  chemical  equilibrium.  GASKET  does,  indeed, 
perform  thia  operation  automatically  by  considering  explicitly  a  fourth  reacr 
tion 


H20  +  CO  J  C02  +  H2 

to  be  in  equilibrium.  This  is  the  so-called  water  gas  shift  reaction. 

In  conclusion,  it  should  be  reiterated  that  the  complexities  in¬ 
volved  in  the  selection  of  base  species  for  mixed  equilibrium  -  nonequilibrium 
systems  are  handled  automatically  by  GASKET  whenever  the  user  selects  all 
program  options  which  utilize  data  built  into  the  program.  In  this  case,  the 
user  needs  to  input  only  the  actual  edge  gas  composition  (see  Section  4). 

However,  if  the  user  chooses  to  input  his  own  kinetic  and  thermochemical  data, 
then  he  roust  carefully  follow  the  procedure  prescribed  in  Section  4  so  that 
GASKET  will  select  the  proper  base  species  and  kir'tic  reactions. 

3.2  FORTRAN  ROUTINES 

In  order  to  permit  an  understanding  of  the  actual  solution  mech¬ 
anics  and  their  relation  to  the  FORTRAN  source  program,  a  brief  description  of 
each  of  the  subroutines  will  be  provided  in  this  section.  For  convenience  the 
routines  have  been  grouped,  where  possible,  under  more  general  headings.  The 
calling  of  these  routines  and  the  principal  iterative  loop  is  controlled  by  the 
mainline  program.  With  the  exception  of  KININ,  SQUEE,  and  ETIMEF  all  routines 
are  called  only  by  the  mainline  program. 

3.2.1  Thermodynamic  Data  Input  Routines 

IMELM:  Reads  elemental  data  information  and  normalizes  ele¬ 

mental  composition  information;  creates  fictitious 
"elements"  when  necessary. 

INPUT:  Reads  species  thermochemical  data  from  DATX  or  input  cards 

selects  base  species,  determines  stoichiometric  coefficients 
of  formation  reactions  of  all  other  species  from  the  base 
species,  flags  condensed  species,  sets  up  first  guesses,  reads 
diffusion  factor  and  fail  temperature  data  (KI. (3)  ^  9)  and 
calculates  molecular  weights. 
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;  ’  I 

BELCH:  Establishes  a  base  species-element  correspondence  table 

which  is  required  if  an  input  element  is  temporily  omit¬ 
ted  from  the  chemical , system  being  analysed. 

DATX:  Contains  built-in  JAN  At'  thermocheniical  data  deck  (see 

Volume  II)  . 


3.2.2  Problem  Setup  and  Initiation 

ZIPIN:  Initializes  and  sets  various  internal  control  flags  and 

reads  shock  data. 

KININ:  Reads  control  card  and  optional  frozen  edge  and  surface 

I  ( 

temperature  arrays;  generates  an  array  of  internal 
control  integers  bo  define  sequence  of  subsequent  cal¬ 
culations  . 

ALPST:  Reads  diffusion  factor  values  if  KRK ( 3 )  =  9;  sets  diffusion 

factors  based  on  exponential  approximation  or  input  values; 
determines  which  elements  are  absent  from  the  system  and 
flags  corresponding  molecular  species  and,  in  particular, 
corresponding  base  species  (see  BELCH) ;  reinitializes 
species  omitted  from  the  prior  solution  and  zeros  those 
to  be  omitted  from  current  solution;  evaluates  conserved 
quantity  parameters  entering  mass  balances. 


3.2.3 


Calculation  of  Errors  and  Error  Derivatives  Used  Within 
Newton-Raphson  Iteration 


THERM:  On  first  iteration  and  after  every  change  of  system  temp- 

|  , 

erature  this  routine  if  called  to  evaluate  molecular  ther¬ 
modynamic  properties  (e.g.,  enthalpy,  entropy,  free  energy, 
specific  heat)  and  the  equilibrium  constants  appropriate 
to  each  formation  reaction.  On  the  first  iteration  cer¬ 
tain  reinitializations  are  performed  including  assigning 
temperature  to  3000°K  if  the  nonconvergent  flag  has  been 
set.  Also  on  first  iteration,  certain  kinetic  rate  data 
are  read  or  generated  by  this  routine  and  certain  key  sum¬ 
mations  are  determined  in  CRECT.  (The  reason  for  these 
functions  inclusion  *n  THERM  relates  to  the  overlaying  of 
subroutines  and  molecular  species  data  in  an  IBM  1130 
version  of  this  program) . 

! 

1 

Initializes  mass  balance  error  equations  and  determines 
contribution  of  bane  specie?  to  these  errors.  Initializes 
matrix  of  partial  derivatives  of  mass  balance  errors  with 


MATL: 
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respect,  to  log  of  base  species  partial  pressures  and  in¬ 
troduces  base  species  contributions  to  these  coefficients. 
Commences  search  for  surface  species.  On  first  iteration 
normalizes  gas  phase  partial  pressures. 

MAT2:  Evaluates  errors  in  formation  reaction  equilibrium  rela¬ 

tion  for  non-base  species.  Introduces  contribution  of 
non-base  species  to  mass  balance  eauations.  For  gases 
the  contribution  of  the  partial  derivatives  of  the  mass 
balance  errors  with  respect  to  the  non-base  species  log 
partial  pressure  is  related  base  species  log  partial 
pressure  via  the  equilibrium  equation  and  introduced  in¬ 
to  the  array  initialized  in  MATl.  For  each  condensed 
species  which  has  been  previously  introduced  into  the 
system  an  additional  equilibrium  relation  is  added  to 
the  mass  balance  equations.  One  additional  condensed 
species  relation  is  accepted  if  equilibrium  indicates  a 
positive  formation  potential  (see  AFMAT  relative  to  this 
species)  Completes  search  for  surface  species. 

MAT3 :  Completes  formulation  of  mass  balance  errors  and  their 

derivatives.  Introduces  condensed  species  coefficients 
into  their  error  equations.  Determines  if  redundancy 
exists  amongst  condensed  species  equilibria.  Removes  any 
redundancy  by  eliminating  either  newly  introduced  species 
(see  MAT2)  or  previously  present  species,  based  on 
equilibrium  errors.  In  latter  case  mass  of  eliminated 
condensed  species  is  distributed  amongst  other  condensed 
species  in  order  to  preserve  mass  balances. 

KINET:  Generates  certain  kinetic  rate  data,  adds  kinetic  rate 

terms  to  mass  balance  equations,  and  evaluates  and  in¬ 
troduces  error  derivatives  into  derivative  matrix.  Re¬ 
arranges  mass  balances  so  that  controlling  reactions  in¬ 
fluence  only  one  mass  balance.  Modifies  coefficients  if 
controlling  reaction  is  approaching  equilibrium  so  as  to 
achieve  linearity  of  equation. 
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3.2.4  Calculation  of  Variable  Corrections 

Upon  inversion  of  the  matrix  of  error  derivatives  and  its  multi¬ 
plication  of  the  errors  (see  RERAY)  a  set  of  corrections  are  obtained  for  temp¬ 
erature,  pressure-molecular  weight  product,  log  partial  pressure  of  the  base 
species  and  the  relative  moles  of  condensed  species.  Before  these  constraints 
can  be  applied,  however,  a  significant  amount  of  messaging  is  required. 

AFMAT:  Immediately  after  the  inversion,  performs  a  series  of 

tests  on  the  corrections  and  if  necessary  imposes  con¬ 
straints  on  equations  and  repeats  inversion.  Examples 
include:  the  newly  introduced  condensed  species  with 
zero  concentration  has  a  calculated  negative  correction 
-  the  corresponding  equilibrium  equation  is  removed  and 
the  correction  set  to  zero;  if  the  temperature  is  cur¬ 
rently  equal  to  a  minimum  or  maximum  value  dictated  by 
phase  change  or  other  discontinuous  phenomena  and  if  a 
contrary  temperature  correction  is  predicted  -  the  temp¬ 
erature  constraint  is  removed  and  the  temperature  frozen; 
if  the  surface  equilibrium  error  and  the  temperature  cor¬ 
rection  are  of  conflicting  signs  -  the  temperature  con¬ 
straint  is  removed  and  temperature  frozen  -  if  an  exces¬ 
sive  negative  change  in  the  pressure  -  molecular  weight 
product  is  predicted  -  the  pressure  constraint  is  removed 
and  pressure  frozen.  Certain  temperature  minima  and  maxima 
are  also  set  to  control  the  convergence  on  surface  equilib¬ 
rium  calculations. 

SCALE:  From  basic  corrections,  corrections  to  log  partial  pressure 

of  non-base  species  are  calculated.  Scans  all  corrections 
and  determines  maximum  damping  factor  which  will  permit 
all  corrections  to  satisfy  certain  constraints.  These 
include  one  order  of  magnitude  increase  and  three  order 
of  magnitude  decrease  in  partial  pressures  of  signifies, it 
species,  and  less  restrictive  constraints  on  less  impor¬ 
tant  species.  The  scale  factor  generated  by  this  routine 
will  subsequently  be  applied  to  all  corrections.  Performs 
correction  of  temperature  and  pressure-molecular  weight 
product. 
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CRECT:  Performs  all  corrections  according  to  scale  factors  cal¬ 

culated  in  SCALE.  Makes  corrections  based  either  on 
linearization  of  mass  balances  (delta  of  mole  fractions) 
or  of  equilibrium  equations  (delta  of  log  of  mole  frac¬ 
tions)  depending  on  relative  importance  in  mass  balances. 
Evaluates  set  of  summations  previously  discussed  rela¬ 
tive  to  THERM  in  Section  3.2.3. 

3.2.5  Calculation  of  Output  Quantities 

PROPS:  Determines  certain  derivative  properties  from  the  in¬ 

verse  matrix  of  error  derivatives  including  the  equilib¬ 
rium  specific  heat  and  the  isentropic  exponent;  calculates 
transport  properties  using  model  based  on  diffusion  fac¬ 
tors.  Calculates  Mach  number  in  shock  stagnation  calcula¬ 
tion,  evaluates  next  pressure  in  iteration  and  assesses 
convergence  at  stagnation  solution;  outputs  calculated 
quantities. 

OUTPT:  Generates  and  outputs  all  additional  terms  required  for 

output  as  displayed  in  Section  5.  Initiates  iteration 
for  shock  stagnation  solution. 

3.2.6  Utility  Routines 

RERAY:  Linear  equation  solver  and  inversion  routine. 

Call  list  is 

N  =  number  of  rows  in  C 
C  =  coefficient  array 
NN  =  number  of  columns  in  C  minus  N 
D  =  set  of  column  vectors  of  length  N 
NNN  =  number  of  column  vectors 
LS  =  column  rearrangement  flag  (zero  implies 
no  rearrangement) 

IS  =  Flag,  -2  at  call  results  in  before  and 
after  display  of  arrays,  <0  at  return 
implies  singularity  encountered  at 
row- IS. 

Returns  with  left  NXN  square  of  C  inverted  and 
this  inverse  multiplied  by  remaining  columns  of 
C  (if  any)  and  D. 
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SWAP:  Replaces  one  array  and  vector  with  another  set,  call 

list  is 

A  array  1  (16  x  16  in  current  listing) 

B  vector  1  (16  in  length  currently) 

SLA  array  2 
SLB  vector  2 

Returns  with  A  and  B  replaced  by  SLA  and  SLB,  respecti¬ 
vely.  SLA  and  SLB  are  unchanged. 


SQUEE:  Prepares  integer  and  alphanumeric  information  from  real 

variables  for  compressed  card  output. 

Call  list  is 

A  real  variable  to  be  converted 
I  first  of  three  storage  locations  for  results 
IDEE  presumed  decimal  location  in  result 
Returns  with  I  (1)  an  integer,  I  (2)  the  sign  of  the 
exponent  (alphanumeric),  I  (3)  the  exponent,  such  that 


A 


Id) 

^IDEE  * 


10 


±1(3) 


ETIMEF 

and 

ETIME:  System  time  routines  available  on  Univac  1108  systems 

tape.  For  other  systems  a  conversion  or  dummy  routine  is 
required.  Call  of  ETIME  sets  zero  time  in  clock.  Sub¬ 
sequent  call  of  ETIMEF  (X)  returns  time  in  seconds  from 
last  call  of  ETIME  as  real  variable  X. 


3.3  DEFINITIONS  OF  SEVERAL  OUTPUT  QUANTITIES 

Several  quantities  appear  on  the  standard  program  output  under 
abbreviated  names.  In  this  section,  these  variables  are  more  accurately  de¬ 
fined.  The  derivative  property  output  3ives  the  following: 

I  3h. 

CP-FROZEN  ■  y  K.  — ,  the  frozen  constant 
i-1  1  wT 

pressure  specific  heat  capacity 
(Btu/lbm°R  or  Cal/gm#K) 
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CP-EQUIL 


the  equilibrium  constant  pressure 
specific  heat  capacity 
(Btu/lbm°R  or  Cal/gm#K) 


DLNM/DLNT  =  §|g2) 


(unitless) 


DLNM/DLNP 


GAMMA 


(unitless) 

§15?)  "  Vcv  '  for  “  ideal  9“ 

's 

(unitless) 


Other  property  output  includes: 

MU1  «  =  £  XjFj ,  a  gas  mixture 

property  which  reduces  to  unity 
for  assumed  equal  diffusion  coef¬ 
ficients  Fj  =  F  =  1.0  (unitless) 

MU2  ■  \>2  a  gas  mixture 

property  which  reduces  to 
the  mixture  molecular 
weight,  ,  for  assumed  equal 
diffusion  coefficient 
Fj  -  F  =  1.0  (grams  of  sys¬ 
tem/mole  of  gas) 

M0L.WT  =  %  ,  the  composite  mixture 

molecular  weight  (see  Equation 
5,  (grams  system/moles  of  gas) 
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HT1L  -  fi  -  T.  Z.h. ,  a  property  of  the 

j-1 

gas  mixture  which  reduces  to  the 
static  enthalpy,  h,  for  assumed 
equal  diffusion  coefficients.  (Btu/lbm) 


I 

CPTIL  =  C  =  V  Z .  — s-=-  ,  a  property  of 

p  i=l  1  dr  p 

the  gas  mixture  which  reduces  to 
the  frozen  specific  heat  capacity 
for  assumed  equal  diffusion 
coefficients.  (Btu/lbm°R  or  Cal/gm°K) 


I 

HTIL*  =  h*  =  £  Zfh.  ,  a  property  of  the 
i«l  1  1 

gas  mixture  which  reduces  to 
HTIL  when  GAMEX  (the  unequal 
diffusion  exponent)  =  1.0,  and 
reduces  tc  static  enthalpy,  h, 
for  assumed  equal  diffusion  co¬ 
efficients.  (Btu/lbm) 


GAMEX  =  y ,  the  weight  factor  exponent  between 

K  and  Z  mass  fractions  used  to  obtain 
the  Z*  values.  (unitless) 

ELEMENTAL  K  MASS  FRACTIONS  =  K^»  usual  definition  of  mass  fraction 

for  element  K  in  the  system.  (gm  of- 
k/gm  system) . 

ELEMENTAL  Z  MASS  FRACTIONS  =  Z£,  elemental  fT  weighted  mass  fraction 

of  element  K  in  the  system,  (gm  of  K/gm- 
system) . 


All  of  the  above  properties  are  calculated  and  output  for  each 
problem  solution  using  the  current  values  of  F^  and  GAMEX,  even  though  the 
particular  option  being  computed  does  not  require  these  valuer .  An  important 
example  is  the  equal  diffusion  surface  equilibrium  option.  Even  though  the 
HTIL*  value  is  computed  and  printed  for  a  GAMEX  value  of  0.667,  the  value  of 
HTIL*  punched  on  output  cards  (for  input  to  a  heat  conduction  program)  will  be 
static  enthalpy,  h,  as  it  should  be. 
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Additional  descriptions  of  the  variables  computed  and  output  by 
the  GASKET  code  are  given  in  References  2  and  3 . 

3.4  PROGRAM  OPERATION  DETAILS 

3.4.1  Program  Storage  Requirements 

The  current  version  of  the  GASKET  program  allows  179  chemical  species 
10  elements,  8  entries  in  the  frozen  edge  temperature  table,  15  entries  in  the 
surface  temperature  table,  and  50  entries  in  the  diffusion  factor  table.  With 
these  limits  approximately  43,000  (decimal)  words  of  storage  are  required.  How¬ 
ever,  roughly  13,600  (decimal)  words  of  this  storage  are  occupied  by  the  thermo¬ 
chemical  data  contained  in  the  routine  DATX . 

3.4.2  Tap:  Requirements 

All  input  and  output  tape  designations  are  set  at  the  beginning 
of  the  main  program  with  one  exception.  The  designations  are 


KIN 

System  input 

(=5) 

KOUT 

System  output 

(=6) 

JAN 

Scratch  (or  save) 

(=18) 

tape,  see  KR(3) 

In  RERAY,  KOUT  is  also  set  to  the  system  output  tape  (=6  currently) . 

3.4.3  Operator  Controls 

No  operator  intervention  control  switches  exist.  All  controls  are 
set  at  the  initiation  of  execution  of  each  problem  with  the  corresponding  con¬ 
trol  card  (see  Section  4). 

3.4.4  Dumps  and  Diagnostics 

Control  column  10  (KRK(10))  controls  the  output  of  requested  dia¬ 
gnostic  data.  The  interpretation  of  this  data  is  not  always  obvious  and  users 
should  consult  with  the  program  authors  if  interpretation  is  mandatory.  This 
same  data  will  automatically  be  obtained  if  more  than  67  iterations  are  re¬ 
quired  for  convergence  and  will  continue  through  the  70th  and  last  iteration. 
Other  diagnostic  complaints  may  be  self-explanatory.  The  failure  to  obtain 
convergence  is  usually  related  to  physical  situations  for  which  no  unique  solu¬ 
tion  exists. 
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3.4.5  Program  Running  Time 

i 

Program  run  time  is  most  conveniently  measured  in  terms  of  the  num¬ 
ber  of  solutions,  and  time  per  solution,  required  for  a  given  problem.  Time 
per  solution  is  conveniently  measured  in  terms  of  the  number  of  iterations,  and 
time  per  iteration,  required  for  a  converged  solution.  Thus,  the  run  time  for 
a  given  problem  may  be  estimated  as: 


Run  Time 


No.  of  Solutions  x 


Iterations 

Solution 


x 


Time 

Iteration 


I 


The  number  of  solutions  depends  upon  the  type  of  problem  being 
solved.  For  example,  when  the  isentropic  expansion  option  is  used  to  calculate 
the  edge  state  and  the  edge  elemental  composition  includes  both  H  and  6,  three 
solutions  are  required  t  obtain  the  final  form  of  the  edge  solution.  Then,  if 
the  built-in  frozen  edge  and  surface  temperature  arrays  are  utilized,  17  addi¬ 
tional  solutions  are  provided,  making  a  total  of  20  solutions  in  this  particular 
problem.  < 


i 


The  number  of  iterations  required  for  a  converged  solution  is,  of 
course,  a  function  of  the  type  of  problem  being  solved.  The  sample  problems, 
Section  5,  should  provide  some  guidance  in  estimating  this  quantity.  Roughly, 
most  solutions,  other  than  the  initial  solution,  converge  in  about  5  to  15  itera 
tions.  1 


The  time  required  per  iteration  is  a  function  of  the  computing 
machine  as  well  as  the  problem  type.  Machines  in  the  speed  range  of  the  Univac 
1108  or  CDC  6600  require  on  the  order  of  10  to  60  milliseconds  per  iteration. 
Again,  the  sample  problems  in  Section  (which  were  run  on  a  Univac  1108)  should 
provide  some  guidance  here.  ! 


I 
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SECTION  4 

I 

INPUT  PREPARATION 

•  This  section  defines  the  format  and  significance  of  each  field  of 

the  input  data  card  deck  for  the  Aerotherm  Graphite  Surface  Kinetics  (GASKET) 
program.  tIn  the  most  general  case,  the  input  consists  of  eight  card  sets.  These 
are: 

I 

1.  CONTROL  CARD 

2.  FROZEN-EDGE  TEMPERATURE  ARRAY 

3.  SURFACE  TEMPERATURE  ARRAY 

4.  ELEMENTAL  COMPOSITION 

5.  DIFFUSION  FACTOR  DATA 

6.  SPECIES  THERMOCHEMICAL  EQUILIBRIUM  DATA 

7.  OBLIQUE  OR  NORMAL  SHOCK  DATA 

8.  REACTION  RATE  DATA 

Card  Set  1  is  a  control  card  which  contains  a  10-element  array, 

i  * 

I^RK(I).  This. array  controls  most  of  the  program  options  and  tells  the  program 

what  to  expect  from  the  remaining  card  sets.  Few  problems  require  all  card 

1  , 

sets  since,  in  most  cases,  the  user  will  wish  to  utilize  the  data  built  into 

I 

the  GASKET  program  to  replace  card  sets  2,  3,  5,  6,  and  8.  The  card  sets  are 
described  in  detail  in  this  section. 


r 

I 

! 
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CARD  SET  1  -  CONTROL  CARD  <1  Card) 

FIELD  1  (Columns  1-10,  Format  1011)  the  variable  array  KRK(10)  which  is 
used  to  control  the  various  program  options. 

COLUMN  1  —  EDGE  THERMODYNAMIC  STATE 


0 

1 


2 

3 


4 


5 


Input 

Static  state  behind  shock,  upstream  conditions 
input  (card  set  7) 

Same  as  1,  except  stagnation  state  behind  shock 
Isentropic  expansion  from  input  chamber  condi¬ 
tions 

Static  state  b>  tind  shock,  upstream  conditions 
calculated  vi,  xsentropic  expansion  from  input 
chamber  conditions  (card  set  7) 

Same  as  4,  except  stagnation  state  behind  shock 


COLUMN  2  --  STATE  OPTION  (pertains  to  Z,  cols.  11-20  of  this  card) 

0  Assigned  temperature 

2  Assigned  enthalpy 

3  Assigned  entropy 


COLUMN  3  --  SPECIES  THERMOCHEMICAL  DATA 


0  9 

0 

1234 

2  4 
34 


9 


Use  built-in  JANAF  data 

Use  built-in  diffusion  factor  data 

Read  from  cards  new  JANAF  data  and,  if  provided, 

diffusion  factor  data 

Print  out  thermochemistry  for  allowable  species 
Store  all  JANAF  data  and  diffusion  factors  read  from 
from  cards  on  tape  (must  be  requested  if  KRK 
(8)  =  4) 

Update  diffusion  factor  data 


COLUMN  4  --  PYROLYTIC  GRAPHITE  SURFACE  SPECIFICATION 


1 

2 

3 


Layer-oriented 

Edge-oriented 

Bulk 


COLUMN  5  —  FROZEN  -  EDGE  TEMPERATURE  ARRAY 


0 

1 

2 


Bypass  this  calculation 
Use  built-in  temperature  array 
Read  the  temperature  array  from  cards  (card  set 
2) 


COLUMN  6  —  SURFACE  TEMPERATURE  ARRAY 


0 

1 

2 


provides  two  equilibrium  surface  calcula¬ 
tions,  for  500°K  and  1000°K,  with  the 
restriction  of  zero  ablation;  this  is  done 
ir.  order  to  provide  tables  for  CMA  and/or 
ASTHMA  which  facilitate  smooth  interpolation. 


Bypass  this  calculation 

Use  built-in  temperature  array 

Read  the  temperature  array  from  cards  (card  set 

3) 

NOTE:  When  KRK(6)>  1,  GASKET  automatically 
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COLUMN  7  --  KINETIC  CONSTANTS 

0  Use  built-in  Arrhenius/Langmuir  model  (see 

Section  2. 2. 2. 2) 

1  3  Read  from  cards  the  pre-exponential  factor, 

B'  (lb  mole  reaction/f t2sec) ,  activation 
energy,  E  (cal/gm  mole),  and  temperature 
exponent,"1* 

23  Read  from  cards  the  partial  pressure  multi¬ 

plicative  factors,  x.  and  K.  (atm“l) 

}m 

COLUMN  8  —  EDGE  GAS  ELEMENTAL  COMPOSITION 

1  e.g.  does  not  contain  H  or  O  (Kinetically 

inert  surface) 

2  e.g.  does  not  contain  H  (reaction  2.  only) 

3  e.g.  does  not  contain  0  (reaction  3.  only) 

4  e.g.  contains  both  H  and  O  (reactions  1., 

2. ,  and  3. ) 

COLUMN  9  —  PUNCHED  CARD  OUTPUT 

0  No  punched  card  output 

1  3  Punch  a  card  appropriate  to  CMA  input  for 

each  convergent  frozen  edge  or  surface 
solution 

23  Punch  a  card  appropriate  to  ASTHMA  input  for 

each  convergent  frozen  edge  or  surface 
solution 


COLUMN  10  —  DIAGNOSTIC  OUTPUT  CONTROL 


0  No  diagnostic  output 

1  Output  a  single  line  of  diagnostic  information 

per  iteration 

j  (Where  j  is  a  digit  greater  than  one)  output 

full  diagnostic  information  for  5 ( j  —  1 )  itera¬ 
tions  and  output  a  single  line  of  diagnostic 
information  for  subsequent  iterations. 


FIELD  2  (Columns  11-20,  Format  F10.5)  Thermodynamic  State  Variable,  Z 


1) 

a) 

If 

KRK ( 2 )  =  0, 

b) 

If 

KRK (2 )  =  2, 

c) 

If 

KRK (2)  *  3, 

2) 

a) 

If 

KRK (1 )  =  0, 

b) 

If 

KRK  (1)  -  1 

c) 

If 

KRK ( 1 )  >  3, 

i)  If  PC 

ii)  If  PC 

’K 


upstream  of  a  shock  wave 
5  pertains  to  the  chamber  thermodynamic  stc 
0,  Z  pertains  to  the  actual  chamber  state 
0,  Z  pertains  to  the  ideal  chamber  state  i 
PC,  Cols.  41-50  of  this  card) 


FIELD  3  (Columns  21-30,  Format  F10.5)  Pressure,  PR  atm 

1)  If  KRK(l)  =  0,  PR  pertains  to  the  edge  thermodynamic  state 

2)  If  KRK(l)  ft  0,  PR  pertains  to  the  actual  chamber  thermodynamic  state 

or  to  the  thermodynamic  state  upstream  of  a  shock  wave 
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FIELD  4  (Columns  31-10,  Format  F10.5)  Edge  or  Shook  Upstream  Thermodynamic 
State,  ZE 

1)  IF  KRK(l)  <  2,  this  field  is  iqnored 

2)  a)  If  ZE  >  0,  it  is  assigned  as  the  pressure,  atm,  at  either  the 

edgo  (KRK  (1)  -  3) ,  or  the  upstream  sido  of  a  shock  wave 
(KRK(l)  -  4  or  5) 

b)  If  ZE  <  0,  it  is  assigned  as  tho  Mach  number  at  either  the 
edgo  (KRK  (1)  -  3)  or  the  upstream  sido  of  a  shock  wave 
(KRK (1)  -  4  or  5) 

FIELD  5  (Columns  41-50,  Tormat  F10.5)  Ideal  Chamber  Pressuro,  PC  atm 

FIELD  6  (Columns  51-60,  Format  F10.5)  Mass  Transfer  Coefficient,  ROUCM 
lbm/ft2coc 

FIELD  7  (Columns  61-70,  Format  F10.5)  Diffusion  Coefficient  Exponent  for 
Unequal  Diffusion  Coefficient  Modol,  GAMER 

1)  If  GAMER  -  0.0,  it  is  automatically  sot  to  2/3  (this  value  is 
appropriate  for  most  calculations) 

2)  If  GAMER  ■  10“5,  tho  diff'*.*ion  model  roduces  to  equal  diffusion 
coefficients 

FIELD  8  (Columns  71-80,  Format  2A4,A2)  Job  Title,  Tile  (3) 


CARD  SET  2  -  FROZEN-EDGE  TEMPERATURE  ARRAY  TFA(I)  (number  of  cards  -  no.  of 
. . .  temperatures  ♦  1) 

1)  If  KRK (5)  •  0  or  1,  skip  this  card  net 

2)  If  KRX(5)  ■  2,  read  in  an  array  c f  no  more  than  eight  temperatures, 

TFA(I)  *K 

CARD ( S )  1,2,...  (Columns  1-10,  Format  F10.5)  one  temperature  per 
card 

FINAL  CARD  blank 

NOTE!  The  TFA(I)  array  built  into  GASKET  contains 
temperatures  from  500  to  4000*K  in  500*K 
increments 


CARD  SET  3  -  SU'tFACE  TEMPERATURE  ARRAY  TSA(I)  (number  of  cards  -  no.  of 
- .  —  ■  - ■■  temperatures  ♦  1) 

1)  If  KRK (6)  -  0  or  1,  skip  this  card  set 

2)  If  l'RK(S)  ■  2,  road  in  an  array  of  no  more  than  fifteen  temperatures, 

TSA(I)  *K 

CARD(S)  1,2....  (Columns  1-10,  Format  F10.5)  one  temperature  per 
card 


FINAL  CARD 


blank 
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NOTE:  The  TSA(I)  array  built  into  GASKET  contains 
temperatures  from  1200  to  4000°K  in  2000°K 
increments.  Also,  whenever  KRK (6)  >  1  two 
surface  equilibrium,  zero-ablation  calcula¬ 
tions  are  performed  for  surface  temperatures 
of  500  and  1000°K. 


CARO  SET  4  -  ELEMENTAL  COMPOSITION  (number  of  cards  -  no.  of  elements  +1) 

CARD  1  (Columns  1-3,  Format  13)  the  number  of  elements  in  the  system 
CARD  2,3,4  .  .  .  (one  card  for  each  element) 

FIELD  1  (Columns  1-3,  Format  13)  the  atomic  number  of  the  element 

FIELD  2  (Columns  4-15,  Format  3A4)  the  name  of  the  element  (for  out¬ 
put  identification  only) 

FIELD  3  (Columns  16-25,  Format  F10.5)  the  atomic  weight  of  the  element 

FIELD  4  (Columns  26-35,  Format  F10.5)  the  relative  amount  of  the 
element  in  the  edge  gas  (this  is  the  same  as  the  relative 
amount  of  the  element  in  the  chamber  or  at  the  upstream  side 
of  a  shock  wave,  since  elements  are  conserved) 

1)  T-isitive  values  are  in  relative  gram- atomic-weights  (or  moles) 

2)  Negative  values  are  in  relative  masses 

FIELD  5  (Columns  46-55,  Format  F10.5)  the  relative  amount  of  the 
element  in  the  surface  material  (the  present  version  of 
GASKET  allows  the  surface  to  be  comprised  of  carbon  only; 
hence,  the  quantity  1.0  must  be  entered  in  this  field  for 
the  element  carbon) 


CARD  SET  5  -  DIFFUSION  FACTOR  DATA  (number  of  cards  =  1/4  x  (7  +  no.  of  data 

— — — —  items,  e.g.,  diffusion  factors,  reference  molecular  weight,  diffusion 

exponent) 

1)  If  KRK (3)  ■  0,  this  card  set  must  be  skipped 

2)  If  KRK(3)  ■  1,2, 3, or  4,  this  card  set  may  be  skipped 

3)  If  KRK (3)  ■  9,  this  card  set  must  not  be  skipped 

CARD  1  (Columns  1-3,  Format  13)  the  total  number  of  data  items  (e.g., 

diffusion  factors,  reference  molecular  weight,  diffusion  exponent) 
to  be  entered 

Diffusion  factors  may  be  specified  for  any  or  all  species 
individually  and/or  diffusion  factors  may  be  calculated  via 


Diffusion  Factor 


Molec . Wt . 
REFM 


iFFA 
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where  REFM  and  FFA  may  be  specified  here.  For  species  for  which 
diffusion  factors  are  not  specified  individually,  diffusion  factors 
will  be  calculated  via  the  above  correlation,  and  if  REFM  and  FFA 
are  not  specified,  diffusion  factors  will  be  calculated  from  the 
above  with  REFM  =23.4  and  FFA  =  0.431.  , 

CARD (S)  2,3,.  .  . 

FIELDS  1,3, 5, 7  (Columns  1-8  ,  21-28  ,  41-48,  61-68  ,  each  Format  2A4) 

The  "name"  of  the  species  for  which  data  is  to  be 
provided,  exactly  as  it  appears  in  columns  73-80  of 
the  first  card  of  the  3  thermochemical  data  cards 
(card  set  6)  for  that  species.  To  specify  values 
of  REFM  or  FFA  in  the  next  field,  enter  the  alpha 
characters  "REFM"  or  "FFA"  respectively  here. 

FIELDS  2, 4, 6, 8  (Columns  9-20,  29-40,  49-60,  69-8Q,  each  Format 
E12.4) 

1)  If  the  name  of  a  chemical  species  was  entered  in 
the  preceding  field,  this  number  is  presumed  to  be 
a  diffusion  factor. 

2)  If  the  name  REFM  or  FFA  was  entered  in  the  preceding 
field,  enter  the  desired  values  of  REFM  or  FFA 
respectively  in  this  field. 

The  diffusion  factors  introduced  into  the  solution  in  this 
manner  replace  either  the  values  set  automatically  (see  card  1  above) 
if  the  species  thermochemical  data  are  read  from  cards  or  tape,  or 
the  values  built  into  the  GASKET  program  (see  Volume  II)  . 

NOTE:  If  KRK (3 )  =  9  and  KRK(l)  =  1,2,4,  or  5,  this  card  set  must 
follow  the  OBLIQUE  OR  NORMAL  SHOCK  DATA  card  set;  if  KRK  (3) 
=1,2,3,  or  4  this  card  set,  if  used,  precedes  the  SPECIES  , 
THERMOCHEMICAL  EQUILIBRIUM  DATA  card  set. 

CARD  SET  6  -  SPECIES  THERMOCHEMICAL  EQUILIBRIUM  DATA  (number  of  cards  =  1  +  3  x 
■  '  no.  of  species) 

1)  If  KRK (3 )  =  0  or  9,  skip  this  card  set 

2)  If  KRK (3)  =  1,2,3,  or  4,  read  in  this  card  set;  there  are  three  of  these 

cards  for  each  molecular,  atomic,  or  condensed 
species;  the  end  of  this  card  set  is  signaled 
by  a  blank  card 


CARDS  1,4,7  .  .  .  Describe  the  elemental  composition  of  the  species  and 
establish  its  name  designation 

FIELDS  1,3,5,.  .  .,  13  (one  for  each  element  in  the  specie  ' 

(Columns  1-3,  7-9,  13-15,.  .  .,  37-j9,  each 
format  F3.0) 

The  number  of  atoms  (of  atomic  number  given 
in  the  following  field)  in  a  molecule  of 
this  species 
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<£E  field  1  is  it-.rc,  this  card  is  presumed 
»:«j  represent  the  end  of  card  set  6) 


FIELDS  2,4,6,.  .  . ,  14.  (cb«  rot  each  «]ec'.»t  in  the  species) 

(Columns  4  -6 ,  10-12 ,  15-12 , .  .  . ,  40-42 , 
each  format  13! 

The.  atomic  numbers  ;«  th%  elements  in  the 
molecule  (tl»e  n-unotu.  of  atoms  of  which 
was  •»ivoj*.  in  the  previous  field) 

( Jf  fluid  i  is  aero,  this  card  is  presumed 
to  L-h  the  first  card  of  card  set  b) 


FIELD  15  (Columns  43-72,  Format  7A4A 1)  the  source  and  date  of  the 

theraochemical  data  for  this  specie.  Used  for  output  only 


FIELD  16  (Columns  73-80,  Format  2A4)  the  name  designation  of  this 
species  (e.g.,  AL203) .  This  variable  id  used  for  output 
and  as  a  means  of  identifying  data  entered  in  card  set  5 


CARDS  2,5,8.  .  .  Lower  temperature  xonge  ♦ hermochemical  data 


FIELDS  1,2, 3, 4, 5  and  6  (Columns  1-54,  Format  6E9.6)  Input  the 

5  constants  (D ,F2,F3,F4,F5,F6)  appropriate 
to  the  J'iwer  temperature  ranqe  of  the  thermo¬ 
dynamic  data  Tor  this  species.  These 
constant?  are  defined  as  follows,  where  T  is 
is  in  *K: 

Fi  »  the  heat  of  formation  of  the  species 
at  299*K  from  the  JANAF  base  state 
(elements  in  most  natural  form  at 
29J»#K)  in  cal/mole 

F2  *  the  enthalpy  change  of  the  species 
from  2S8°K  to  3000®K  in  cal/mole 

F3.F4  and  F5  arc  defined  by  a  cur fit  to 
the  heat  capacity  at  constant  pre.-^uia  of 
the  form: 


f,  »  F3  i  F4  T  *  in  cal/mole  "K 

*• 

F5  *=  the  entropy  of  the  species  at  3000°K 
in  cal/mole  *K 


FIELD  7  (Columns  61-56,  Format  F6.0)  the  upper  limit  of  the  lower 
temperature  range  in  °X 


FIELD  8  (Column  67,  Fermat  II)  the  phase  spec i f ication : 
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1  signifies  gaseous  species 

2  signifies  solid  species 

3  signifies  liquid  species 

The  only  phase  combination  allowed  in  one  three  card  set  is  solid- 
liquid  in  Which  case  a  2  and  3  would  appear  in  fields  8  of  cards  2 
and  3,  respectively. 


CARDS  3,6,9.  .  .  Upper  temperature  range  thermochemical  data 

These  cards  are  the  same  as  cards  2,5,8.  .  .  except  use 
constants  for  the  upper  temperature  range  and  field  7  is 
ignored 


LAST  CARD  OF  CARD  SET  6  MUST  BE  BLANK 

The  end  of  the  species  thermochemical  data  is  signified  by  a  blank  card. 
Hence,  the  last  card  of  card  set  6  must  be  a  blank  card. 


See  Section  3.1.3  for  additional  comments  on  the  ordering  of  the  species 
thermochemical  data  and  the  selection  of  base  species. 

NOTE;  If  this  data  is  read  in  on  cards,  the  species  at  the  top  of 
the  deck  must  be  in  the  following  orders  C0? ,  H~,  H_0,  CO, 

C* ,  C-H-,  HC1 ,  HF,  ....  Furthermore,  the  atomic  number  of 
C*  must  be  106,  rather  than  the  actual  value  of  6.  Finally, 
all  co  densed  phase  species  (except  C*)  must  follow  the  gas 
phase  tpecies. 


CARD  SET  7  -  OBLIQUE  OR  NORMAL  SHOCK  DATA  (1  card) 

1)  If  KRK(l)  =  0  or  3 ,  skip  this  card  set 

2)  If  KRK(l)  =  1,2,4,  or  5,  this  card  set  must  be  read  in 

FIELD  1  (Column  1,  Format  II)  Input  the  variable  KVH  which  determines 
how  the  velocity  ahead  of  the  shock,  the  static  enthalpy  ahead 
of  the  shock,  and  the  total  enthalpy  are  calculated  (since 
total  enthalpy  =  static  enthalpy  +  1/2  velocity2) 

If  KVH  =  1,  velocity  is  calculated  from  specified  total  and 
static  enthalpy 

If  KVH  =  2,  total  enthalpy  is  calculated  from  specified  static 
enthalpy  and  velocity 

If  KVH  =  3,  static  enthalpy  is  calculated  from  specified  total 
enthalpy  and  velocity 


FIELD  2  (Columns  2-10,  Format  F9.4)  The  velocity  variable,  UR  (this 
field  is  ignored  if  KVH  =  1) 
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If  UR  /  -9999,  then  UR  is  assigned  to  SVL,  the  velocity  upstream 
of  the  shock  in  ft/sec 

If  UR  =  -9999,  the  velocity  calculated  in  the  prior  solution 
is  assigned  to  SVL 


FIELD  3  (Columns  11-20,  Format  F10.4)  The  total  enthalpy  variable, 

HTR  (this  field  is  ignored  if  KVH  =  2) 

If  HTR  /  0,  then  HTR  is  assigned  to  HCH,  the  total  enthalpy 
of  the  system  in  BTU/lbm 
If  HTR  ■  0,  HCH  is  unchanged 

FIELD  4  (Columns  21-30,  Format  F10.4)  the  static  enthalpy  variable, 

HSR  (this  field  is  ignored  if  KVH  =  3) 

If  HSR  /  -9999,  then  HSR  is  assigned  to  SHI,  the  static  enthalpy 
upstream  of  the  shock  in  BTU/lbm 

If  HSR  =  -9999.,  the  static  enthalpy  calculated  in  the  prior 
solution  is  assigned  to  SHI 


FIELD  5  (Columns  41-50,  Format  F10.4)  the  shock  angle  (0.0  for  a 
normal  shock)  in  degrees 


CARD  SET  8  -  REACTION  RATE  DATA  (number  of  cards  =  5  x  MT  at  maximum) 

1)  If  KRK (7 )  =  0,  skip  this  card  set 

2)  If  KRK (7)  =  1  or  3,  the  first  subset  of  this  card  set  must  be  read  in 

3)  If  KRK (7)  =  2  or  3,  the  second  subset  of  this  card  set  must  be  read  in 

First  Subset:  MT  cards 

CARDS  1,2,  .  .  .  ,  MT 

FIELD  1  (Columns  1-10,  Format  E10.4),  TKF (M)  the  pre-exponential 
factor  for  the  Mth  reaction,  lb  mole/ft2sec  atm 

FIELD  2  (Columns  11-20,  Format  E10.4),  EAK(M),  the  activation  energy 
for  the  Mth  forward  reaction,  cal/gm  mole 

FIELD  3  (Columns  21-30,  Format  E10.4),  EXK(M) ,  the  temperature  exponent 
for  the  Mth  reaction 

Second  Subset:  4  x  MT  cards  (for  maximum  of  10  base  species) 

CARDS  MT  +  1,  MT  +  2;  MT  +  5,  MT  +  6;.  .  .  5  x  MT  -  3,  5  x  MT  -  2; 

FIELDS  1-IS  (Columns  1-10  x  IS,  Format  8F10.6,  will  require  2  cards 

if  IS  >  8)  inhibiting  specie  partial  pressure  coefficient, 
PSI  (J,M) ,  J  =  1,  IS,  atm-1  (see  Section  2. 2. 2. 2) 
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CARDo  MT  +  3,  MT  +  4;  MT  +  7,  MT  +  8;  .  .  .  5  x  MT  -  1,  5  x  MT 

FIELDS  1-IS  (Columns  1-10  x  IS,  Format  8F10.6,  will  require  2  cards 

if  IS  >  8)  inhibiting  specie  partial  pressure  coefficient 
FKK(J,M),  J  =  1,  IS,  atm-1  (see  Section  2.2. 2.2) 

NOTE:  MT  =  no.  of  reactions 


1) 

IF 

KRK (8)  =  1, 

MT  =  0,  kinetics  ignored 

2) 

If 

KRK (8)  =  2, 

MT  =  1 

M  =  1 

C*  +  C02  *■  2CO 

3) 

If 

KRK (8)  =  3, 

MT  =  1 

M  =  1 

2c*  +  h2  *  c2h2 

4) 

If 

KRK  (8)  =4, 

MT  =  4 

M  =  1 

c*  +  h2o  J  co  +  h2 

M  *  2 

C*  +  C02  J  2CO 

M  =  3 

2C*  +  H2  £  C2H2 

M  =  4 

H20  +  CO  *  C02  +  H2 

IS 

=  no.  of  base  species 

1) 

If 

KRK (8)  =  1, 

kinetics  ignored 

2) 

If 

KRK (8)  =  2, 

J  =  i  co2 

J  =  2  CO 

J  =  3  C* 


J  >_  4  remaining  base  species 

3)  If  KRK (8 )  =3,  J  =  1  H2 

J  =  2  C* 

J  =  3  C2H2 

J  ^  4  remaining  base  species 

4)  If  KRK (8)  =4,  J  =  1  C02 

J  =  2  h2 

J  =  3  H20 

J  =  4  CO 

J  =  5  C* 

J  =  6  HC1  or  HF 

J  _>7  remaining  base  species 


TO  TERMINATE  A  JOB 


KRK(l)  =  8  on  final  card 
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MULTIPLE  JOBS 

1)  Several  jobs  may  be  loaded  back  to  back,  with  a  control  card 
signifying  the  start  of  each  new  job: 

CONTROL  CARD  #1 

ELEMENTAL  COMPOSITION  #1 

additional  data  as  required  for  job  #1 

CONTROL  CARD  #2 

ELEMENTAL  COMPOSITION  #  2 

additional  data  as  required  for  job  #2 

KRK ( 1 )  =  8 


2)  If  jobs  following  the  first  job  have  the  same  elemental  composi¬ 
tion  as  the  first  job,  Card  Set  #4  (ELEMENTAL  COMPOSITION)  may  be 
replaced  by  a  single  blank  card  for  all  jobs  after  the  first: 

CONTROL  CARD  #1 

ELEMENTAL  COMPOSITION  #1 

additional  data  as  required  for  job  #1 

CONTROL  CARD  #2 

blank  card 

additional  data  as  required  for  job  #2 
KRK (1)  =  8 


c 
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SECTION  5 
SAMPLE  PROBLEMS 

Presented  in  this  section  are  five  sample  problems  which  were  run  on  a 
Univac  1108  digital  computer.  An  attempt  was  made  to  utilize  almost  all  of  the 
options  available  in  the  GASKET  program.  For  each  sample  problem,  the  following 
is  presented: 

•  A  brief  description  of  the  nature  of  the  problem  and  solution 

•  A  listing  of  the  input  data  deck 

•  A  few  typical  pages  of  the  program  output 
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Sample  Problem  1 

In  this  problem,  ablation  rates  at  the  throat  of  a  rocket  nozzle  are  cal¬ 
culated.  The  propellant  gases  are  comprised  of  seven  elements,  including  both 
hydrogen  and  oxygen.  The  program  options  which  utilize  built-in  physical  data 
are  selected;  hence,  only  nine  cards  of  input  data  are  required  for  the  complete 
sequence  of  calculations.  The  ideal  chamber  conditions  (propellant  flame  tempera¬ 
ture  of  3764 °K,  ideal  pressure  of  68.07  atm)  are  input.  After  redefining  the 
elemental  set  to  include  the  fictitious  element  "GRAPHITE",  subsequent  calculations 
performed  by  the  program  include  determination  of  the  actual  chamber  state  (54.42  atm) 
and  an  isentropic  expansion  to  the  throat  (unity  Mach  no.).  The  boundary  layer 
edge  state  at  the  throat  is  then  recomputed  with  a  redfined  elemental  set  which 
includes  the  fictitious  element  "HOT  DRY  ICE";  this  is  done  for  the  purpose  of 
the  subsequent  surface  kinetics  calculations .  The  frozen  edge  state  is  computed 
for  the  built-in  array  of  eight  temperatures,  and  the  surface  ablation  rate  is 
computed  for  the  built-in  array  of  fifteen  surface  temperatures.  Also,  a  punched 
card,  with  format  appropriate  for  CMA  input,  is  output  for  each  of  the  frozen 
edge  and  surface  solutions.  Molecules  whose  concentrations  are  smaller  than 
machine  capacity  appear  with  zero  mole  fractions  in  the  output. 


3001110*10  176*.  b*.«2  -1.0  68. 02  0. 8*9  0.  PR  o.  I  WELATIVt  ELEMENTAL  COMPOSITIONS*  ATOMIC  uls/UNIT  MASS 

7  AT. NO.  ELEMENT  ATOMIC  »T  EQuE  OAS  SURFACE 
1  ht  IMOGEN  1.008  3.3182  I  HYDROGEN  1.00600  .0330707  .0000000 
*  BERYLLIUM  8.012  l.*26S  *  BEPYLLIUM  <9.01200  .01*2128  .0000000 
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Sample  Problem  2 

The  shockwave  option  of  GASKET  is  utilized  in  this  problem.  The  boundary 
layer  edge  state  is  specified  as  the  stagnation  condition  behind  a  normal  shock 
wave  with  an  upstream  Mach  number  of  two.  The  shock  is  assumed  to  exist  in  the 
divergent  section  of  a  nozzle  and,  thus,  conditions  upstream  of  the  shock  are 
calculated  via  an  isentropic  expansion  from  input  ideal  chamber  conditions.  Since 
the  edge  gases  contain  both  hydrogen  and  oxygen,  all  three  kinetically-controlled 
surface  reactions  are  possible,  requiring  that  the  edge  state  be  recalculated 
in  terms  of  an  elemental  set  which  contains  the  fictitious  element  "HOT  DRY  ICE". 
Two  temperatures  are  input  for  both  the  frozen-edge  and  surface  temperature  arrays. 
Finally,  the  diffusion  factors  are  updated  for  sixteen  selected  species. 
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Sample  Problem  3 

In  this  problem,  ablation  rates  at  tho  throat  of  a  rocket  nozzle  are  again 
calculated.  The  actual  chamber  state  is  input,  from  which  an  isentropic  expansion 
to  the  nozzle  throat  is  performed  in  order  to  obtain  the  boundary  layer  edge  state. 
Since  the  propellant  gases  do  not  contain  oxygen,  only  one  of  the  three  possible 
kineticelly-controlled  surface  reactions  can  occur  (reaction  3.,  Section  2. 2. 2. 2) 
and,  therefore,  the  fictitious  element  need  not  be  introduced  in  the  edge  state 
definition.  No  frozen  edge  solutions  are  calculated,  and  an  array  of  three  sur¬ 
face  temperatures  is  input. 


GtMPHlTC  SU«f»CE  KINETICS  (WSKETI  SXuTION  PWN.  J  C-KpmITC  So»f«C(  KINETICS  P*OP.  ]  )H«IH'M 


r-i 


U 

91 

X 


* 


I 

wt  n 
W  • 
I  w 

— 

§ 

»- 

5 


3 

a 


o 

o 


5 


5 

a 

® 


ix 

o 


i 


o 


l/> 

1/1 

V 

a 


o 

* 


«j 


o  -» 

o  « 

o  o 


N 


O  O 
«/» 

4  O 

1  s 

N  O 
o  If 
u  • 
,/l  rw 
/  V 

UJ  *P 


a  a 

*  * 
4  UJ 


»  —  U 

GOO 
•  ♦  |  © 
O  n  o  «n 

lA  O  M  <J‘ 

•  >  M  4 

* /i  n<Mi  • 

4 

o  •  •  • 


X  oooc 

u.  o  o  o« 

o  o  o  c 

«J  O  O  O  € 


ON  U.  4 

a  £  *  «  x  £ 
iflU uuuuu 


a  o  <oo  o  o  o 

O  '*>00000 


3  0  — 


tf  * 


u  u.  u. 
o  4  X 

a*  unmz  cn 
7i  u  Xooooo 


£ 

T 

®  .X  «■» 

J(®H^ 
S  Jl4  I 

©  ►*  ©Uj 
T.  -J  M 

I 

<  I  1 

a  *■  ►>  _j 

J  Qt  •-  5 
4  O-iU 
X  a:  J\  U4 
*-  •-  2  _J 
*  .2  U>  C 
WWOl 


•  <  * 

a  400 
N.  4  9  OO 

4>  noo 
wo  "OO 
.**•—00 


>  O  Cf 

>  O  K 

>o« 


21S 

UK*N 

i/nv^n 

p.  r  — 


!l 


* 

3 


•3  or  ♦  ►  - 

waMa« 

6  my  IM 


4  -  « 

w  O  w 

•  •I  — 
W  04)  40 
k  a)«  49 

—  y  <7-  t,n 

e/1  <n  o  e/i  • 

o  —  c  —« 


</»  O  o 

v/>  y  OO  O 

W  IP  O  O  o 

2  7  0  0  o 

a,  woo  • 


y>  £ 


O  O  O  O  O  O 

o  o  o  o  o  o 
o  o  o  o  o  o 
o  o  o  o  o  o 

3  0  0  0  0  0 


ON  U.  # 

uj  x  Mr 

x  n  *  n  r  n 
l/»  o  u  u  u  u  u 


•  ry  1/1 

3  0  0 

►  I  • 

©  —  « 

I  tj  JO 
—  —  /> 


*? 

*5  . 

\  fl  J4  S-* 

I"  XSZz 

g  •  -  •  • 

u. 

•  *  - 

i?  5  :  r 

j!  5HH 

—  O  I  I 

8  IMM 

L;  1  i  : 

ij?  i*H 

8C  i 


UJ  o  o 

!  Vi 


ni 


o  •  * 

“  4“  4-  P- 

4  •  •  i 

4-  P- 

41  2 

>- 

*  i 

o  o  o 

•/too  — 

2  0  0  0  — 

O  o  O  3 

—  o  o  3 

►  •  •  O 

O  i/I 

«  9 

rt  a 

In.  —  —  3 

o  o  — 

in  j  »  -x 

<n  9  o  m  x 

«  « •*  - 

X  (M  IT*  l/t  3 

ON  2  - 

-  -.t  §  2 

fp-  U4 
(A  4 

<  u-  a  * 

O  O  X  UJ  UJ 

/  O  O  —  —  — 

I  I  P-  —  E/1 

J  Otfl  > 

*  *•  *  Z  A 

►-  —o  o  o 

2  fM  n  —  O 

u  N  4  u  UJ 

g  •  •  3  ^ 

3  £  3 

UJ  1/1  u 


®  O  4ff 

e/i  /m  9  i/i  • 

4  N  2  —  « 

O  •  •  • 


i/i  & 

« 

O  IT 

UJ 

•A  C 
yi 

;  ? 
N  O 
Q  A 

UJ  • 
</>  — 
2  *■ 


a 

UJ 

o 


3V 

t-  3  X  — 
TH  tU 
3  J* 
I 

t  •  a 
>•  o 
a  ►  ►  3 
j  a.  p-  5 
<m  o  —  U 
IB  /U 
p-  ►-  2  _J 
Z  2  WO 
uj  uj  O  X 


a  oxooooo 
u.  0-00000 
o  M  O  o  O  O  O 
3  O  — OOOOO 
0—00003 
. . .  • 


u  u.  u 
yj  2  t 

£•  XM 

AOIUOOOU 


•  I  •  I 

X  OXOOOOIM 
u  11-000041 
o\oooo  4 
UJ  o  —  O  O  O  O  e/i 
J  «  — ooo  o  > 


uj  w  n 

a  N  M  4  X 

/I 4UUUDI 


Sample  Problem  3,  Sample  Output,  Continued 


r 


5-55 


MO 

Si  —  ~ 

r  <x 


•  •  l 

x  yj  a*  o  «  n 

*-  H>  n  O  4  K 

<  **  -C  o  o  f> 

i/tiNOO  • 
9  O  ^  O  A 

I  !“* 


.  900^00 
ooo^oo 

jOOO  «oo 


I  > 

4  n 

5  - 


X  *X  .J  o 

it 

*a  **bC 


$  3 


3  J  O 

O  —  ♦ 

2  z 
-  o  e 


-j  Z  o 
>  uj  • 

-  >S«  T 

-  O  X 

«  a  /> 

»  I*,  o 

-  i  — 

a  a  • 

«*J  u 


•  X  T>  _i  X 
J  «  —  —  -3 
y  ®  ^ 

75^  r  ^ 


i  /i  •-  « 
I  3 

X  Z  *  •  »* 

w  O  —  _/  c 

U  A*  C  4 

•  X  ♦ 


go  —  fM  4> 

»/>  •©  ?  — 

-*n  t  o 

U.  >  AJ  — 

7  •  * 

3  4  © 

loo 
•  I 
*•  o  o 

»-  a  o«4 

1103® 


F: 


p  o  o 

I  ii 

7  I  N 
©  «  » 
j  -«  7  "• 
►.  4>  4 


o  o 
S<  O  o 

2  UOO 

5  o  o 


1  (V>  « 
7  <X 
•x  o  > 


o  o 

t  0  3 

J  ©  i 

4  O  <7 

K  ~  »  I 

Z  «X 
y  *X  * 


y»  « 

1 2 


>-o- 

-O  1 
4  *-» 

a  • 


*-  —  i 


CS 


a  a 


¥  ^ 1 


X  ~ 
JO  X 

—  —  o 

st» 

JS3 

uuu 
w  c 
4  a  r 
jji 

A*b 

~  Z  -J 
*  * 

jj  i 


y.  *A 
ItOO 
3  4  X 
J^xit 


o 

uJ  o  o 
7  o  o 
zoo 


©  o  o 

iii* 

o  o  ■>  o 

vA  fx  4'  O  • 
4  flifi 

J  t  »  I 


1  " 
d  a 

o 

*8 


-i  a 
»-  x 

<  ^ 


X 

®  X  n 

v  J  y.  C 
3  x  n  «3 

>-  D  X  — 

■tin 

*  •  t 


Z  Z  .y 
uJ  U  3  X 


*  ’ 

X  • 

LO 

4  3 

9  3! 

-J  Z 

ui  *  u* 

0>-  9 

4  O  7 


131 

o  1 

◦  uj 

*■*  * 
*-  ■  3 
Z  Ui 

uj  *  /1 


3  *  Z 

4  4  o 

o  o 


OX  U  i 

OX  At  I 

a  -X  #^I  X 

/»  O  o  o  O  o  o 


O  f\j 

•  •  I 

X  OK  31000 
U.  04)04000 
040  OOO 
JO  «0  A  ©  ©O 
J  O  40  1  O  30 


U  k  4. 

J  7  X 

1  «  J  *»(X  1  Xfti 

Zoloouuu 


a  0KO9PO9 

U.  A  C  O  T  ©  o  » 
4>  4'  ©  n  it  o  — 

y  X  mo  41  4  00 

jta  rfO  mmO  O 


p  a 

Jik  W 


Xi  U  U. 

*?5 


J  c  -  ® 
>A  "X  *-  4 
a*  X  tZ 


•  I  4 

**  *  o  ^  K 
-00*4 
yiho«  • 
3  NO  *  « 

i ... 


1000700 

4000  too 

O  O  O  4)  o  o 
y  o  o  o  o  o  o 
j  o  O  O  KlO  O 


o  —  —  . 


~  4 

:  x  ©  —  45 

ja>i- 
1  '•»  J  N 


4)  O 
>  Zn  o 
OJ'O 


zoo 

zoo 

UJ  o  o 
3  0  0 

z  •  • 


tj  IX  4  J 

y  I  1  *x  r 

1  ix  tux 

A  J  U  U*  u>  O  V 


a  s> 
Z  I 
-J  — 

o  o 


►*  u>  ( 

7  il 


*  s  * 

in'- 


UJ  z  o 
»  UJ  ♦ 
MN  t 

►  00 

«  a  *x 

»l»  X 

tl  • 

JJ  u 


y  I  • 

O  — 

•  1  f  JX 


3  -x  m  •» 

3  Z  O  •  f» 

f  O  4)  -J  4> 

O  f  O  4> 


»-  -O  -X 
3  O  i/>  »X  C 

o/»  3- 

—  »  f  o 

T 


a  c  —  4 
x  x  o  d  a 

ui  UJ  O  X  — 

a  *-  >  4 


o  o 

1  o  o 

COO 


u.  —  -• 

o  o 
«Z  I  ft 

J»  4  N. 

4  x>  o 

X  <X9  «c 

o  «x 

>4  0  7 


I  I 

»  o  <7 

4>  7 

■A  fX 
1  -X  4 


•  o 
o 

•  r 


?32 

►  o  • 

a  a 

3  *- 


X  u. 

o  U  UJ 

.  ^3 

y 

<  a  t 

j  y  t 

fl)l  J 
~  Z  -J 
4  4 
J  X  ■ 
u  ►- 
4  » 

u.  <A 
1  'A  o 
o  4  r 

X  X 


00c 
•  •  •  4 
Z  C7  C  *1 
i7  ►» 
on  ©  4> 
u'  ^  OP  —  • 
4  <X  t  z 
3  •  •  • 


*j  y 

r  o 

o 
O  O 


Z  O 
y  o 
U  7 


x  x  n 

-J  I  *-  *- 
X  JU  t 

3  \  \  J 

*-  3  I  •— 

s>-  ry 

X  J  * 

I 

I  I  X 

>■  4 

a  »  >-  -1 

j  1  k  5 

4  O  -•  %J 

1  a  i/'y 
*-  »-  Z  -1 
74410 

y  y  1.'.  > 


1  O  K  o  0000 
»  o c 00000 
O  4  0  0  0  0  0 
jO  MOOOOO 
J0^305*»3 


U  1 
O  uJ 
►*  -J  M 

<3 

*-  X  o 

y  h  / 

o  z 

«  l/l  UJ 
3  Ui  o 

3  71  z 

4  4  3 
I  o  O 

uJ 

•-  1  I 
4  O  J 
►*  U.  u. 


•  I  • 


—  <x 


1  J]  KOU1P04I 
u.43  O  o  O  —  O  ►• 
o  <■  o  o  S>  O  T> 
y«  xOOOu  4 
JK  •*  O  AX  O  X 

o  . . 


O  ly  u> 

—  _J  _ 

105 

iilt 

r - 


I 


Sample  Problem  3,  Sample  Output,  Continued 
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Sfflpl.  Problem  J.  Semple  Output,  Concluded 
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Sample  Problem  4 

Again  ablation  rates  at  the  throat  of  a  rocket  nozzle  are  calculated. 

The  edge  state  is  computed  via  an  isentropic  expansion  from  input  actual  chamber 
conditions.  Since  the  propellant  gases  do  not  contain  hydrogen,  only  one  of  the 
three  possible  kinetically-controlled  surface  reactions  can  occur  (reaction  2., 
Section  2. 2. 2. 2)  and,  therefore,  the  fictitious  element  need  not  be  introduced 
in  the  edge  state  definition.  Unlike  the  four  other  sample  problems  presented 
in  this  section,  the  equal-diffusion  coefficient  mode  is  utilized  here.  No 
frozen  edge  solutions  are  computed,  and  an  array  of  five  surface  temperatures  is 
input. 
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Sample  Problem  4,  Sample  Output,  Continued 
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Sample  Problem  5 

In  this  problem,  all  options  which  require  the  input  of  thermophyaical 
properties  are  utilized.  The  edge  thermodynamic  state  (2005. 81*K  and  12.778  atm) 
ia  input  directly.  Diffusion  factors  for  sixteen  selected  species  are  input, 
along  with  a  complete  JANAF-type  thermochemical  data  deck.  Pre-exponential  factors; 
activation  energies,  temperature  exponents,  and  inhibiting  species  partial  pres¬ 
sure  coefficients  are  input  for  each  of  the  four  possible  kinetically-controlled 
reactions  (the  three  reactions  in  Section  2. 2. 2. 2  and  the  water  gas-shift 
reaction  mentioned  in  Section  3.1.3).  Finally,  two  frozen-edge  temperatures 
and  two  surface  temperatures  are  input  to  be  used  in  lieu  of  the  corresponding 
built-in  temperature  arrays. 
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Sample  Problem  5,  Sample  Output,  Concluded 
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